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ABSTRACT 


Omega  is  a  worldwide  VLF  navigation  system  which  utilizes  phase 
differencing  techniques.   The  VLF  signals  of  Omega  are  perturbed  by 
phase  uncertainties  caused  by  factors  such  as  diurnal  effects  which 
are  quite  predictable  and  other  factors  which  are  not  so  predictable. 
This  dissertation  models  the  effective  phase  velocity  uncertainties 
which  result  from  the  less  predictable  phase  uncertainties  as  jointly 
Gaussian  perturbations  to  the  phase  velocities.   This  approach  of 
modeling  VLF  propagation  for  Omega  has  not  previously  been  taken.   The 
more  predictable  factors  are  incorporated  into  determination  of  the 
mean  values  of  the  phase  velocities.   This  model  allows  any  interfre- 
quency  and  interstation  correlation  to  the  phase  disturbances  of  Omega 
to  be  readily  taken  advantage  of  to  improve  the  estimation  performance. 
From  this  model,  the  maximum  likelihood  position  estimation  equations 
are  derived  for  Omega  with  an  arbitrary  number  of  frequencies.   Monte 
Carlo  simulation  results  are  reported  for  four-frequency  Omega  as  SNR 
and  phase  velocity  uncertainty  are  varied.   These  results  allow  quanti- 
tative error  prediction  for  Omega  utilizing  different  fourth  frequencies 
to  be  made  for  the  first  time.   Several  fourth  frequencies  which  can 
be  easily  implemented  by  the  Omega  stations,  to  place  less  stringent 
requirements  on  ambiguity  resolution  techniques,  are  examined.   The 
requirement  for  resolving  the  inherent  lane  ambiguity  of  Omega  arises 
when  application  of  Omega  to  search  and  rescue  systems  such  as  the 
planned  worldwide  Global  Rescue  Alarm  Net  (GRAN)  is  considered.   Four- 
frequency  experimental  data  was  processed  with  the  maximum  likelihood 
estimator  and  the  results,  which  demonstrate  the  potential  of  the 
estimator,  are  presented. 
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I.   INTRODUCTION 

The  Omega  navigation  system  is  being  developed  with  the 
intent  of  achieving  a  worldwide  navigational  aid.   It  is 
anticipated  [1],  that  the  eight  station  Omega  system  will  be 
fully  operational  by  1976.   There  are,  at  this  writing,  two 
fully  operational  stations  radiating  the  prescribed  10  kw 
of  power  at  the  three  present  Omega  frequencies;  10.2  kHz, 
11  1/3  kHz,  and  13.6  kHz.   These  stations  are  located  in 
North  Dakota  and  Norway.   The  exact  locations  are  contained 
in  Table  I.   The  Trinidad  station,  which  will  eventually  be 
replaced  by  the  station  in  Liberia,  is  operational  at  re- 
duced power;  and  the  Hawaii  station  is  expected  to  be  fully 
operational  at  full  power  by  1975. 

The  choice  and  number  of  frequencies  for  the  Omega  sys- 
tem has  been  predicated  on  the  use  of  Omega  as  a  navigation 
aid  by  ships  and  aircraft  which  would  have  dead  reckoning 
(DR)  and  alternate  fix  capabilities.   In  the  recent  past, 
however,  the  United  States  Navy  [2]  and  the  United  States 
Air  Force  [3],  along  with  several  other  organizations,  have 
been  investigating  the  use  of  Omega  for  applications  such 
as  search  and  rescue  (SAR) .   The  requirements  of  a  SAR  Sys- 
tem, using  Omega  as  a  position  location  information  source, 
are  much  more  demanding  than  those  of  navigation  systems. 
In  general,  no  a  priori  position  information  is  available 
at  a  SAR  incident  location.   On  the  basis  of  what  information 
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TABLE  I 
OMEGA  STATION  LOCATIONS 


Station Lat  i  tude Longi  tude 

Aldra,  Norway  66°25*15"  N  13°09'10"    E 

Trinidad**  10°42'06.2"  N  61°38'20.3"  W 

** 
Liberia   (proposed) 

Haiku,  Hawaii  21024*16.9"  N  157 ° 4 9 ' 52 . 7 "  W 

La  Moure,  North  Dakota  46°21,57.2"  N  98o20*08.8"  W 

Reunion  20°58'26.5"  S  55°17'24.2"  E 

Trelew,  Argentina  43o03'12.5"  S  65°11'27.7"  W 

Australia  (proposed) 

Tsushima,  Japan  34°36'53.3"  N  129°27'12.5"  E 

Forestport,  New  York*  43°26'40.9"  N  75°05'09.8"  W 


Information  for  this  table  was  obtained  from  the 
Defense  Mapping  Agency,  Washington,  D.  C. 

** 

The  Trinidad  station  is  expected  to  be  replaced  by 

a  station  in  Liberia. 

*** 

The  Forestport  station  remains  only  as  an  experi- 
mental station  and  is  not  part  of  the  Omega  system. 
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can  be  received  at  the  SAR  incident  and  relayed  to  rescue 
forces,  the  SAR  position  must  be  determined. 

A.   INTRODUCTORY  BACKGROUND 

The  first  proposal  to  use  satellite  relayed  Omega  infor- 
mation for  SAR  purposes  was  by  Samek  and  Pike  [4],   The 
principle  was  taken  from  the  Omega  position  locating  experi- 
ment (OPLE)  proposed  by  NASA's  Goddard  Space  Flight  Center. 
The  OPLE  was  carried  out  in  1967[5]  utilizing  the  ATS-1  and  3 
satellites  and  demonstrated  that  Omega  navigation  signals 
could  be  received  at  a  remote  site  and  retransmitted  via 
satellite  to  a  ground  center  which  uses  these  signals  to 
determine  the  location  of  the  retransmitter.   The  success  of 
this  experiment  engendered  the  advanced  OPLE  concept.   Ad- 
vanced OPLE  was  to  use  a  UHF  satellite  link  frequency  (402  MHz) 
instead  of  the  VHF  link  frequency  (149  MHz)  of  OPLE,  smaller 
transceivers  (platforms),  random  access  techniques  for  the 
satellite  link,  and  an  ambiguity  resolution  technique  for  the 
Omega  information.   Omega  ambiguities  are  described  later  in 
this  chapter.   These  changes  to  OPLE  allowed  for  much  more 
diversified  anticipated  usage  of  the  concept.   The  random 
access  technique  (many  users  using  a  single  satellite  at 
arbitrary  times)  obviated  the  requirement  to  interrogate  the 
platforms  and  thus  they  could  be  much  smaller  and  lighter. 
This  allowed  for  planned  application  of  the  concept  to  tasks 
such  as  tracking  of  animals  for  zoological  purposes.   Reso- 
lution capability  for  the  ambiguity  problem  would  mean  no 
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a  priori  position  information  would  be  required  to  locate 
the  transceiver.   The  advanced  OPLE  concept  was  not  carried 
out  because  of  fund  limitations  and  the  interest  of  Naval 
Air  Test  Center  (NAVAIRTESTCEN) ,  Patuxent  River,  Maryland 
in  using  the  equipment  for  their  GRAN  program. 
1.   GRAN 

In  1969,  NAVAIRTESTCEN,  Patuxent  River,  Maryland, 
proposed  the  Global  Rescue  Alarm  Net  (GRAN)  as  a  worldwide 
SAR  system  based  on  the  advanced  OPLE  concept.   GRAN  will 
utilize  hand-held  transceivers  called  SAR  communicators 
(SARCOM's)  to  receive  Omega  signals  and  relay  them,  shifted 
in  frequency  and  occupying  reduced  spectrum, via  SAR  satel- 
lites (SARSAT's)  to  ground  SAR  centers  (SARCEN's),  Figure  1. 
At  the  SARCOM,  unit  identification  information  is  added  to 
the  spectrum.   In  1970,  NAVAIRTESTCEN  carried  out  experi- 
ments utilizing  low  transmitted  powers  and  demonstrated  the 
feasibility  of  low  power  SARCOM's  for  GRAN  [6].   The  100  kHz 
band  of  frequencies  from  406.0  MHz  to  406.1  MHz  was  assigned 
for  use  in  a  SAR  satellite  link  such  as  GRAN  at  the  World 
Administrative  Radio  Conference  for  Space  Telecommunications 
(WARC-ST)  in  Geneva  in  1971  [2]. 

It  is  anticipated  by  the  originators  that  GRAN  will 
serve  all  countries.   GRAN  will  be  available  to  any  military 
and  civilian  users  who  purchase  a  SARCOM.   In  keeping  with 
this  concept,  the  cost  of  the  SARCOM  units  must  be  kept  to 
the  minimum  possible  level.   This  constraint  prohibits  pre- 
processing of  the  received  Omega  information  in  the  SARCOM. 


12  (   x 


55 
W 
CJ 

si 

CO 


& 


55 
O 

M 

a 

to 


5 

CO 


to 


Q 


§ 

4-1 
CO 
>. 

CO 


55 
O 
H 

a 

CO 

< 


<J 


H 
M 


s 

O 

I 

CO 


o 


0) 

n 

3 
60 


13 


\ 


In  order  to  meet  the  portability  requirements  necessary  so 
pleasure  craft  and  perhaps  even  hikers  can  carry  the  units, 
weight  and  size  of  the  SARCOM  units  must  also  be  kept  to  a 
minimum.   Studies  by  NASA  [7]  have  demonstrated  the  feasi- 
bility of  the  following  characteristics  for  the  SARCOM's: 

Volume  1000  in 

Weight  1  kg 

Output  Power  5  watts  maximum 

Frequency  406.0  -  406.1  MHz 

Bandwidth  2.5  kHz 

Transmission  Period    3  min.  broadcast,  repeated  as 

battery  power  allows 

Identification  Data    36  (bits)  -  sufficient  for 

nine  digit  number  such  as 
a  social  security  number 

The  SARCOM  units  will  operate  sequentially  in  three 
modes  when  activated.   The  modes  are: 

Mode  1.  Acquisition 

Mode  2.  Identification  Data 

Mode  3.  Omega  Data 

In  mode  one,  the  SARCOM  will  transmit  only  an  internally 
generated  Acquisition/Recognition  (A/R)  tone  in  the  allotted 
band  of  frequencies  for  a  period  of  about  12  seconds.   The 
maximum  power  capability  of  the  unit  will  be  concentrated 
in  this  tone  to  allow  the  SARCEN  phase-locked-loop  receivers 
to  acquire  the  retransmitted  tone  and  switch  to  a  narrow 
tracking  bandwidth  for  signal-to-noise  ratio  improvement. 
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The  second  mode  is  then  commenced  automatically  and  SARCOM 
unit  identification  information  is  transmitted  digitally  by 
phase-shift-keying  the  A/R  tone  with  a  code  preasslgned  to 
the  SARCOM  unit.   The  second  mode  will  consist  of  about  six 
seconds  of  transmissions.   In  the  final  mode,  the  received 
Omega  information  is  retransmitted  along  with  the  A/R  tone. 
The  total  transmitted  power  of  the  SARCOM  will  be  divided 
between  the  A/R  tone  (approximately  .5  w)  and  the  retrans- 
mitted Omega  signals  (remaining  power  equally  divided). 
This  third  mode  is  expected  to  last  for  approximately  180 
seconds,  which  would  allow  for  18  of  the  ten  second  Omega 
format  periods  (see  the  following  section)  to  be  relayed. 
A  possible  spectrum  of  signals  to  be  retransmitted  if  four- 
frequency  Omega  is  adopted  and  the  fourth  signal  is  at  10.88 
kHz  is  shown  in  Figure  2.   The  signals  in  the  spectrum  cor- 
respond to  frequency  shifted  versions  of  the  Omega  signals; 
signal  number  one  corresponds  to  the  10.2  kHz  signal,  signal 
two  to  the  10.88  kHz  signal,  signal  three  to  the  11  1/3  kHz 
signal  and  signal  four  to  the  13.6  kHz  signal.   The  actual 
number  of  different  frequency  Omega  signals  and  the  options 
available  in  choosing  these  frequencies  is  discussed  in  de- 
tail later  in  this  thesis.   The  center  frequency  of  the 
spectrum  in  Figure  2  is  referred  to  only  as  /  .   This  is  done 
since  to  utilize  the  entire  100  kHz  bandwidth  of  the  SARSAT 
as  well  as  to  avoid  concurrent  identical  spectra  to  be  trans- 
mitted from  several  SARCOM's,  the  2.5  kHz  SARCOM  spectra 
will  have  to  be  spread  throughout  the  406.0  MHz  -  406. 1  MHz 
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band   of  the  SARSAT.   This  may  be  accomplished  by  assigning 
different  frequency  A/R  tones  to  different  users.   Reference 
[8]  contains  additional  discussion  in  this  area.   The  spac- 
ing of  the  retransmitted  Omega  signals  relative  to  the  A/R 
tone  will  be  accurately  controlled  by  the  SARCOM.   Thus, 
when  the  phase-locked-loop  receiver  in  the  SARCEN  locks  on 
and  tracks  the  A/R  tone,  the  Omega  signals  may  be  recovered 
with  essentially  all  of  the  frequency  drift  and  doppler 
shifts  which  affect  the  spectrum,  after  its  creation  in  the 
SARCOM,  eliminated. 

The  SARSAT  system  to  be  utilized  for  GRAN  will  con- 
sist of  a  system  of  three  to  four  geosynchronous  satellites 
and  additional  polar  orbiting  satellites.  The  exact  number 
and  types  of  orbits  of  the  satellites  is  still  being  con- 
sidered by  NAVAIRTESTCEN.  A  discussion  of  the  requirements 
for  the  satellite  system  as  well  as  link  calculations  is 
contained  in  [8], 

The  SARCEN  system  will  consist  of  several  SARCEN's 
in  order  to  allow  some  SARCEN  to  continuously  be  in  the 
coverage  of  each  SARSAT.   Each  SARCEN  will  search  the  pre- 
assigned  band  of  frequencies  with  a  set  of  phase-locked-loop 
(PLL)  receivers.   The  bandwidth  of  the  PLL  in  the  search 
mode  will  be  larger  than  when  tracking  a  SARCOM  relayed 
message  in  order  for  the  acquisition  time  to  be  small. 
When  an  A/R  tone  has  been  acquired,  the  bandwidth  of  the  PLL 
will  be  decreased  to  improve  the  SNR  from  the  PLL  and  allow 
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lock  to  be  maintained  after  the  power  of  the  A/R  tone  is 
decreased  for  the  third  mode  of  SARCOM  transmission. 

After  the  SARCEN  has  acquired  a  retransmitted  SAR 
message,  identified  the  sending  SARCOM  from  the  identifica- 
tion data,  and  recorded  the  retransmitted  Omega  signals, 
the  position  of  the  SAR  incident  must  be  determined.   The 
optimum  method  for  locating  the  retransmitting  SARCOM  posi- 
tion is  determined  by  this  thesis. 
2 .   Omega 

a.   History  of  Omega 

Following  World  War  II,  the  MIT  Radiation  Labora- 
tory which  had  been  actively  developing  military  navigation 
aids,  was  closed  [9].   Loran  A  and  Loran  C,  as  presently 
used,  were  results  of  their  work.   A  small  group  from  Divi- 
sion 11  of  the  Laboratory  moved  to  Cruft  Laboratory,  Harvard 
University.   In  1947,  Professor  J.  A.  Pierce  of  this  group 
proposed  Radux  [4]  as  a  navigation  system.   The  technique 
of  Radux  was  very  similar  to  Loran  except  that  the  trans- 
mitted signals  were  phase  modulated  with  a  200  Hz  signal. 
Because  of  the  extremely  stable  phase  characteristics  of 
signals  in  the  proposed  frequency  range  of  40  kHz  to  50  kHz, 
it  was  expected  that  the  modulation  could  be  used  in  phase 
comparison  with  quite  accurate  results. 

The  results  of  experiments  with  Radux  by  the 
Navy  Electronics  Laboratory  indicated  that  the  accuracy  of 
Radux  could  resolve  ambiguities  which  would  result  from  a 
10  kHz  signal.   This,  along  with  the  fact  that  VLF 
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frequencies  were  relatively  unusable  for  communications  pur- 
poses, resulted  in  a  composite  system  Radux-Omega  which 
radiated  the  40  kHz  signals  of  Radux  interrupted  by  10  kHz 
bursts  of  a  coherent  signal.   The  10  kHz  signal  yielded 
nominally  eight  nmi  (nautical  miles)  ambiguities  which  could 
be  resolved  by  the  200  Hz  modulation  of  the  40  kHz  signal. 

It  was  soon  found  that  this  system  was  not  as 
good  as  could  be  obtained.   The  10  kHz  signal  had  a  much 
longer  usable  range  than  the  40  kHz  signal  and  could  thus 
have  much  longer  baselines  (the  great  circle  connecting  two 
Omega  stations)  thus  improving  the  accuracy  of  the  system 
markedly . 

b.   Present  Omega 

Several  different  frequencies  in  the  10-14  kHz 
range  were  used  in  experiments  with  Omega.   The  frequencies 
decided  upon  for  final  implementation  were  10.2  kHz,  11  1/3 
kHz,  and  13.6  kHz  common  to  all  Omega  stations  and  other 
frequencies  unique  to  the  stations  to  be  transmitted  for 
other  than  phase  difference  navigation.   A  total  of  eight 
stations  were  decided  upon  to  achieve  worldwide  coverage 
with  a  station  redundancy  sufficient  to  assure  an  adequate 
fix  anywhere  in  the  world.   The  three  common  frequencies  are 
radiated  as  on-off-keyed,  t ime- division- mult iplexed  signals 
in  a  ten  second  format  as  shown  in  Figure  3.   The  nominally 
one-second  pulses  are  all  separated  by  0.2  seconds.   Each 
station  transmits  the  common  frequencies  in  the  same  sequence, 
but  starting  at  a  different  time. 
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Omega,  being  a  phase  differencing  navigation  sys- 
tem is  necessarily  a  hyperbolic  system.   That  is,  the  isophase 
contours  of  two  stations  (lines  on  the  earth  for  which  the 
phase  difference  of  two  common  frequency  signals  from  two 
stations  is  constant)  form  a  set  of  hyperbolas  shown  in 
Figure  4,  which  has  as  foci  the  two  stations  radiating  the 
signals.   As  shown  in  Figure  4,  two  adjacent  hyperbolic  lines, 
henceforth  called  "lines  of  position"  or  LOP,  separated  by  a 
distance  d   on  the  baseline  can  be  shown  to  be  spaced  by 
d   =  d-  esc  (0/2)  off  the  baseline  [10].   The  angle  6  is  that 
angle  formed  by  the  two  great  circles  passing  through  the 
stations  and  a  point  midway  between  the  adjacent  LOP  off  the 
baseline.   The  spacing  of  two  LOP  off  the  baseline  of  two 
stations  then  will  be  equal  to  the  spacing  of  these  LOP  on 
the  baseline  multiplied  by  esc  (8/2). 

Because  it  is  a  phase  difference  system,  Omega 
has  naturally  occurring  ambiguities.   On  the  baseline  of 
two  stations,  the  same  phase  difference  of  a  common  frequency 
f  radiated  from  the  stations  will  be  seen  every  -r-  A  where  X 
is  the  wavelength  of  the  frequency  f.   This  results  in  lanes 
between  unambiguous  LOP  of  width  —  X  on  the  baseline.   Off 
the  baseline  this  spacing  will  be  -y  X    esc  (0/2). 

This  concept  is  easily  extended  to  multiple  fre- 
quency transmissions  by  considering  a  signal  of  frequency  f 
which  is  a  common  clock  frequency  to  all  Omega  stations  and 
is  utilized  to  create  the  K  transmitted  signals  f.  (i  =  1,  2, 
.  .  .,  K) .   In  the  present  Omega  system  K  is  three.   In  order 
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for  phase  comparison  of  the  common  frequencies  from  the  sta- 
tions to  be  useful  in  locating  a  remote  receiver  site,  the 
frequency  f   at  all  stations  are  maintained  very  nearly 
phase  coherent.   The  ith  Omega  signal  to  be  transmitted  is 

created  from  division  of  f   by  an  integer  n,.   The  baseline 

r       j  &     i 

lane  width  associated  with  frequency  f.  will  be  -r-X    =  v  /2f  , 

1        i  2  s    s    s 

where  f   is  the  largest  common  multiple  frequency  of  the 
transmitted  frequencies.   This  is  the  smallest  frequency 
which  can  be  obtained  by  frequency  differencing  techniques 
from  the  K  transmitted  frequency.   The  value  of  v  ,  the 
effective  phase  velocity  of  f  ,  is  determined  by  the  phase 
velocities  of  the  K  signals.   However,  v   is  nominally  equal 
to  c,  the  velocity  of  light  in  vacuum,  so  we  can  get  a  nomi- 
nal value  for  77X   by  assuming  v   ~  c .   The  term  "unambiguous 
2  s  s 

lane"  will  be  used  in  this  thesis  to  refer  to  -r-X    .   The 

2  s 

nominal  baseline  size  of  the  unambiguous  lane  may  be  deter- 
mined by  the  following  procedure.   The  integer  divisors, 

K)  are  first  separated  into  the  products 

of  powers  of  their  prime  number  components  (see  example  which 
follows).   The  product  of  the  largest  powers  of  all  prime 
numbers  which  occur,  then  gives  the  divisor  n   which  yields 


n±  (i  =  1,  2, 


f   from  f   =  f  ,/n  .   Baseline  unambiguous  lane  width  is  then 
s        s     f   s 

~  c/2f   =  c  •  n  /2f  .   The  present  value  of  f  (816  kHz)  used 
s         s    r  r 

by  the  Omega  stations  allows  a  simple  rule  of  thumb  to  find 

a  nominal  baseline  unambiguous  lane  width  easily  from  n  . 

s 


-=-A  (nmi)  ~ 
2  s 


161,875.2(nmi/sec)«  n, 
2  •  816,000  Hz 


=.0992  n  ~  n  /10   (1) 
s    s 
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The  present  three  frequencies  transmitted  by  the  Omega  sys- 
tem; 10.2  kHz,  11  1/3  kHz,  and  13.6  kHz  yield  n  =  2 "  •  3 2  •  5 
=  720  since  n1    =  2  * •  5  =  80,  n2  =  23-32  =  72,  and  n3  = 
22  •  3  *5  =  60.   Thus  from  (1)  the  baseline  width  of  the 
unambiguous  lanes  is  nominally  72  nmi. 
3 .   VLF  Propagation 

Since  the  Omega  system  depends  upon  phase  comparison, 
unexpected  phase  shifts  in  the  propagation  path  cause  posi- 
tion estimation  errors.   The  causes  of  the  anomalous  phase 
shifts  must  therefore  be  considered.   For  the  Omega  system, 
frequencies  in  the  10-14  kHz  region  have  been  chosen  because 
of  the  extremely  long  ranges  which  VLF  signals  propagate  and 
because  of  the  remarkable  phase  stability  of  these  signals. 
The  propagation  problem  at  VLF  has  been  studied  both  experi- 
mentally and  theoretically  by  many  researchers  for  many  years. 
References  [11,  12,  13,  and  14]   contain  very  detailed  and 
thorough  discussions  of  VLF  propagation. 

In  general,  for  distances  less  than  several  Megameters 
(Mm)  from  the  transmitter,  the  easiest  theoretical  analysis 
of  VLF  propagation  can  be  obtained  by  a  combination  of  ground 
and  sky  wave  contributions.   The  ground  wave  provides  the 
field  which  would  be  seen  in  the  absence  of  the  ionosphere 
and  its  attenuation  with  distance  from  the  transmitter  is 
generally  attributable  to  earth  curvature  [11].   The  sky  wave 
is  the  transmitted  wave  which  is  reflected  in  the  earth  iono- 
sphere cavity.   The  contribution  to  the  total  field  from 
skywaves  can  be  determined  by  summing  all  of  the  skywaves 
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present  at  any  distance  from  the  transmitter  using  ray 
theory  [11],   At  distances  beyond  a  few  Mm,  the  number  of 
skywave  contributions  becomes  unwieldy  and  waveguide  theory 
of  propagation  becomes  much  easier  to  use. 

Waveguide  theory  for  VLF  is  generally  used  for  mak- 
ing theoretical  predictions  for  Omega.   Propagation  predic- 
tions (to  be  discussed  later)  do  attempt  to  predict  phases 
of  signals  at  ranges  where  ground  wave  effects  are  signifi- 
cant; however,  Omega  is  considered  most  predictable  and 
thus  most  reliable  at  distances  in  excess  of  2-3  Mm  from 
the  transmitters.   At  these  distances,  single  mode  waveguide 
theory  seems  to  describe  the  propagation  quite  well. 

Waveguide  theory  of  propagation  describes  the  earth 
and  ionosphere  as  forming  a  spherical  waveguide.   The  height 
of  the  ionosphere  above  the  earth,  nominally  90  km  at  night 
and  70  km  during  the  day,  is  on  the  order  of  several  wave- 
lengths for  the  10-14  kHz  signals  of  Omega  and  thus  wave- 
guide theory  would  seem  to  apply  if  the  boundaries  of  the 
waveguide  could  be  adequately  described.   It  is  generally 
accepted  [15],  that  the  propagation  effects  causing  phase 
variation,  which  is  what  is  of  concern  for  a  phase  comparison 
system  such  as  Omega,  can  be  divided  into  two  categories, 
temporal  and  spatial.   Among  the  temporal  variations,  the 
diurnal  variations  have  the  most  dramatic  affect  and  are 
probably  one  of  the  most  predictable  factors.   This  variation 
is  in  general  accounted  for  by  using  different  effective 
phase  velocities  for  each  frequency  for  dark  and  for  lighted 
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transmission  paths.   Professor  J.  A.  Pierce  of  Harvard  Uni- 
versity has  conjectured  that  adequate  position  estimation 
using  Omega  may  be  attainable  using  the  single  waveguide 
mode  phase  velocities  given  in  Table  II.   These  values  were 
obtained  from  [16]  and  from  personal  communication  with 
Professor  Pierce.   Other  temporal  variations  are  long  term 
dependence  on  solar  cycle,  lunar,  and  annual  variations. 

The  spatial  variations  of  importance  appear  to  vary 
with  propagation  bearing,  earth  conductivity  over  the  trans- 
mission path,  geomagnetic  latitude,  and  the  ionosphere  con- 
ditions.  The  magnetic  field  has  an  apparent  marked  effect 
on  the  ionosphere  and  thus  creates  phase  variations  depen- 
dent upon  the  transmission  path.   These  spatial  variation 
causes  are  not  independent  quantities  and  thus  must  be 
treated  together  [15]. 

The  Defense  Mapping  Agency  attempts  to  account  for 
all  of  the  above  phase  variation  sources  as  well  as  other 
lesser  known  sources  in  the  Omega  Propagation  Correction 
Tables  which  are  published  regularly  for  the  10.2  kHz  and 
13.6  kHz  Omega  signals.   The  computer  programs  which  accom- 
plish this  are  described  in  [15].   Additionally,  force 
fitting  of  the  corrections  to  actual  data  is  continually 
carried  out.   The  use  of  these  tables  is  simple;  knowing  the 
date  and  the  approximate  receiver  location,  phase  correction 
factors  to  be  added  to  the  observed  phase  of  a  known  fre- 
quency signal  from  a  known  station  are  looked  up  in  pre- 
published  tables.   After  these  corrections  are  made,  the 
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TABLE  II 
PREDICTED  FIRST  MODE  PHASE  VELOCITY  RATIOS 


f (kHz) 


v/c  Day 


v/c  Night 


10.2 
11  1/3 
13.6 
10.88 


1.00270 

1.00139 

.99965 

1.00187 


.99960 
.99869 
.99751 
.99902 
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phase  difference  for  a  particular  station  pair  and  frequency 
choice  is  computed  and  the  corresponding  LOP  is  determined 
from  specially  prepared  Omega  charts.   Figure  5  shows  a 
sample  section  of  an  Omega  chart.   These  charts  are  prepared 
based  upon  a  phase  velocity  1.00261  times  the  speed  of  light 
in  vacuum  and  thus  the  phase  corrections  attempt  to  correct 
the  phase  to  this  effective  phase  velocity  in  determining 
the  receiver  location. 

B.   PROBLEM  STATEMENT  AND  APPROACH 

This  dissertation  determines  the  optimum  procedure  to 
establish  an  Omega  receiver's  location.   In  order  to  accom- 
plish this  task,  a  model  for  the  Omega  problem  was  first 
formulated  in  which  the  effective  phase  velocities  of  the 
transmitted  Omega  signals  were  assumed  known.   Since  true 
phase  velocities  of  the  signals  are  continuously  varying 
along  the  transmission  path  as  ionospheric  densities  and 
many  other  factors  change,  effective  phase  velocity  is  used. 
This  velocity  is  that  which  satisfies  the  equation 


v  = 


* 


(2) 


where  d  is  the  distance  measured  on  the  globe  which  the 
signal  travels,  to  is  the  angular  frequency  of  the  signal, 
and  (f>  is  the  observed  phase  in  radians  of  the  signal  at  a 
distance  d  trom  the  transmitter. 

Utilizing  this  assumption  and  the  assumptions  about  the 
noise  which  will  be  given  in  Chapter  II,  the  maximum  likeli- 
hood LOP  estimation  equation  for  two-station,  multiple 
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frequency  Omega  was  derived.   This  derivation  is  presented 
in  Appendix  C.   Using  this  result  and  the  effective  phase 
velocities  for  single  mode  propagation  given  in  Table  II, 
some  f our- f requency  experimental  Omega  data  was  processed. 
The  data  was  collected  preliminary  to  the  four  frequency 
Omega  experiment  described  in  detail  in  Chapter  VI.   The  re- 
sults of  this  effort  were  no  lane  errors  (defined  in  Chapter 
IV)  for  transmission  paths  in  daylight,  but  a  significant 
number  of  lane  errors  for  dark  or  partially  dark  transmis- 
sion paths.   The  conclusion  which  was  reached  then  was  that 
the  Table  II  effective  phase  velocities  were  not  sufficient 
predictions  of  phase  velocities  actually  being  experienced. 
Accounting  for  diurnal  effects  only,  in  predicting  effective 
phase  velocities,  and  assuming  single  mode  propagation  was 
not  apparently  justified. 

A  further  search  into  the  literature  revealed  that  at 
reasonably  long  ranges  from  VLF  transmitters,  predicted 

phase  velocities  using  propagation  prediction  methods  dis- 

-4 
agreed  with  experimental  results  by  a  few  parts  in  10   [17], 

This  indicated  that  the  best  available  methods  of  predicting 
the  disturbances  to  the  phase  of  VLF  signals  for  typical 
distances  could  yield  errors  in  phase  differences  of  signals 
from  two  stations  on  the  order  of  one  radian  or  more.   This 
realization  led  to  the  model  which  was  used  in  this  disser- 
tation . 

Since  phase  is  related  to  effective  phase  velocity  by 
(2),  phase  disturbances  can  equivalently  be  considered  as 
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due  to  effective  phase  velocity  uncertainties.   The  many 
effects  which  cause  phase  uncertainties  and  thus  effective 
phase  velocity  uncertainties  imply,  as  a  result  of  the  cen- 
tral limit  theorem,  a  model  of  the  phase  velocities  of 
jointly  Gaussian  random  variables.   This  new  approach  not 
previously  considered  in  the  literature  allows  any  correla- 
tion between  the  phase  disturbances  of  signals  from  the 
same  station  or  different  stations  to  be  taken  into  account 
in  the  estimation  process.   The  effective  phase  velocity 
can  be  predicted  as  accurately  as  necessary  or  possible  for 
the  different  frequency  signals  and  these  values  used  for 
the  mean  values  of  the  jointly  Gaussian  phase  velocities. 
The  entries  in  the  covariance  matrix  of  the  phase  velocities 
reflect  the  uncertainties  of  these  predictions  and  the 
correlation  of  the  uncertainties.   This  is  the  first  model 
for  position  estimation  which  allows  advantage  to  be  taken 
of  the  correlation  between  different  frequency  signals  from 
an  Omega  station. 

References  [18,  19,  and  20]  reveal  a  significant  correla- 
tion between  phase  perturbations  of  Omega  signals  transmitted 
from  a  single  station.   The  interstation  correlation  is  in 
general  considered  to  be  negligible  [19].   If  some  appro- 
priate standard  deviations  of  the  phase  velocities  on  the 

-4 
order  of  a  few  parts  in  10    are  assumed  and  the  correlation 

between  all  signals  are  assumed  known,  the  probability  den- 
sity functions  of  the  assumed  jointly  Gaussian  phase  velo- 
cities are  completely  known.   These  parameters  must  be 
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experimentally  determined.   Professor  J.  A.  Pierce  of  Harvard 
University  is  presently  investigating  the  correlation  of 
phase  disturbances  of  Omega  signals. 

Utilizing  this  new  approach  and  the  assumption  of  addi- 
tive white  noise,  which  is  justified  in  Chapter  II,  the  maxi- 
mum likelihood  estimation  equations  for  Omega  were  derived 
(Chapter  IV)  for  high  SNR.   The  performance  of  the  resultant 
estimator  for  various  fourth  frequencies  which  could  be 
easily  added  to  the  present  three- frequency  Omega  format  was 
investigated  by  Monte  Carlo  simulation.   The  results  of  these 
simulations  (see  Chapter  V)  allow  for  the  first  time,  a 
quantitative  measure,  in  terms  of  expected  errors,  of  the 
performance  of  Omega  with  different  fourth  frequencies  which 
could  be  added  to  Omega  to  yield  larger  unambiguous  lanes. 
The  maximum  likelihood  estimator  was  used  to  process  experi- 
mental four- frequency  data  (Chapter  VI)  and  the  performance 
was  found  to  conform  in  general  to  the  simulation  results. 
Preliminary  comparison  of  the  estimator  performance  indicates 
it  is  superior  in  terms  of  lane  errors  (defined  in  Chapter 
IV),  the  most  important  consideration  for  SAR  systems  parti- 
cularly, to  any  other  known  estimation  procedure. 

The  equation  for  the  maximum  likelihood  estimator  re- 
quires prior  estimates  of  the  doppler  shifts,  signal  ampli- 
tudes and  noise  spectral  density  heights  for  the  Omega 
signals  received  by  the  SARCOM.   The  maximum  likelihood  esti- 
mation equations  for  these  quantities  are  contained  in 
Appendices  A  and  B.   The  results  of  application  of  the  signal 
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amplitude  and  noise  density  estimation  equations  to  experi- 
mental data  are  presented  in  Appendix  B  also.   These  equa- 
tions allow  optimum  estimates  of  signal-to-noise  ratios  to 
be  found  for  Omega  signals.   The  optimum  method  for  determin- 
ing SNR  of  Omega  signals  has  not  been  presented  before. 

A  four  frequency  experiment  utilizing  a  satellite  link 
is  planned  by  Naval  Air  Test  Center,  Patuxent  River,  Maryland 
for  late  1974.   A  prototype  portable  transceiver  is  being 
constructed  and  will  be  transported  to  various  sites  in 
Europe  and  the  eastern  United  States  for  relay  of  four  fre- 
quency Omega  signals.   The  LES-6  satellite  will  be  used. 
The  processor  for  the  ground  station  which  has  been  selected 
for  use  is  the  maximum  likelihood  estimation  program  pre- 
sented in  Appendix  E  of  this  dissertation. 

C.   ORGANIZATION  OF  THESIS 

Chapter  II  of  this  thesis  presents  the  mathematical 
model  of  the  estimation  problem.   The  assumptions  necessary 
to  describe  the  problem  so  that  a  solution  can  be  obtained 
and  their  justification  are  presented. 

Chapter  III  addresses  the  resolution  of  the  ambiguity 
problem  which  arises  in  a  differencing  system.   The  size  of 
the  unambiguous  lanes,  determined  by  the  number  of  signals 
and  the  choice  of  their  frequencies,  determines  the  require- 
ments from  an  ambiguity  resolution  technique.   Several 
potential  techniques  and  their  anticipated  capabilities, 
which  will  thus  determine  frequency  choice  for  a  SAR  system, 
are  presented. 
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Chapter  IV  contains  the  derivation  of  the  maximum  likeli- 
hood estimate  of  the  LOP  resulting  from  two  stations  each 
transmitting  an  arbitrary  number,  K,  of  common  frequencies. 
Since  there  are  only  eight  locations  in  the  Omega  format  for 
the  nominally  one  second  pulses  of  different  frequencies,  K 
will  necessarily  be  less  than  or  equal  to  eight  unless  more 
than  one  frequency  is  allowed  to  occupy  a  pulse  by  time  shar- 
ing.  The  resultant  estimation  equation  when  K  is  set  equal 
to  four,  the  number  of  frequencies  presently  being  considered 
by.  the  managers  of  GRAN,  is  presented.   The  estimation  equa- 
tion which  resulted  when  deterministic,  known  phase  velocities 
were  assumed  is  addressed  in  Appendix  C. 

Chapter  V  presents  the  results  of  Monte  Carlo  simulation 
which  was  carried  out  for  several  easily  implemented  fourth 
frequency  choices,  an  expected  range  of  variances  of  the 
phase  velocities,  and  a  range  of  SNR.   The  simulations  were 
carried  out  to  aid  in  selection  of  a  fourth  frequency  for 
Omega  which  will  allow  an  ambiguity  resolution  technique  to 
resolve  the  unambiguous  lane  problem.   The  trade  off  in  per- 
formance in  terms  of  estimation  error  which  is  experienced 
as  unambiguous  lane  size  is  increased  as  a  function  of  the 
necessary  parameters  is  presented.   Although  the  addition  of 
one  frequency  to  the  Omega  format  is  expected  to  be  accept- 
able to  the  program  managers  as  well  as  the  Omega  station 
host  countries,  additional  frequencies  and/or  modulation  of 
a  frequency  has  been  stated  to  be  unacceptable  by  the  Omega 
program  office  and  is,  therefore,  not  considered. 
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Chapter  VI  presents  the  results  of  processing  four-fre- 
quency Omega  data  from  two  stations  with  the  maximum  likeli- 
hood estimator.   The  data  for  the  experiment  was  collected 
by  Texas  Instruments  Corporation  of  Dallas,  Texas  for 
NAVAIRTESTCEN  in  the  fall  of  1973.   The  fourth  frequency 
which  was  transmitted  by  the  two  experimental  Omega  stations 
(Trinidad  and  Forestport,  New  York)  in  addition  to  the 
standard  three  Omega  signal  frequencies  was  10.88  kHz.   A 
discussion  of  the  results  and  an  explanation  of  the  errors 
which  were  experienced  is  presented.   The  estimation  was 
done  first  by  assuming  the  effective  phase  velocities  given 
in  Table  II  for  single  mode  propagation  and  second  by  using 
propagation  correction  factors  provided  by  the  Defense  Map- 
ping Agency.   A  comparison  of  the  results  is  presented. 

Chapter  VII  presents  the  derivation  of  the  multiple  sta- 
tion estimation  equation  and  the  prospects  of  being  able  to 
implement  its  computer  solution.   An  alternative  approach 
to  the  multiple  station  problem  which  requires  less  computa- 
tion time  but  is  nonoptimum  is  presented. 

Chapter  VIII  presents  the  conclusions  drawn  from  the  pre- 
vious chapters.   Several  appendices  contain  important  methods 
and  results  which  are  not  required  in  the  general  flow  of 
the  chapters. 


35 


1   \ 


II.   MATHEMATICAL  MODEL 

Consider  two  Omega  stations  each  transmitting  K  frequen- 
cies occupying  K  of  the  eight  nominally  one  second  positions 
in  the  predetermined  format  shown  in  Figure  3.  Let  Y(t)  be 
a  2K  dimensional  column  vector,  where  the  first  K  components 
represent  the  signals  received  from  station  one  and  the  last 
K  components  represent  the  signals  received  from  station  two 
We  shall  represent  Y(t)  as 


Y(t)  =  S(t)  +  N(t) 


(3) 


where  S(t)  represents  the  signal  portion  of  the  received  vec- 
tor and  N(t)  represents  an  additive  noise  process.   We  can 
express  the  ith  component  of  S(t)  in  (3)  as 


s.(t)  =  A.  cos  [(to.  +  A.)t  -  4>  .  +  9,]     (4) 
l        x         x     x      Ti     i 


where  A.  is  the  amplitude  of  the  ith  signal,  u) .  is  the  angu- 
lar frequency  of  the  ith  signal,  A.  is  the  doppler  frequency 
shift  of  the  ith  signal,  $ .     is  the  propagation  phase  delay 

of  the  ith  signal  and  0.  is  the  contribution  of  other  than 
b  x 

propagation  delay  to  the  phase  of  the  ith  signal.  The  com- 
ponents of  the  2K  dimensional  signal  vector  will  be  ordered 
so  that  the  angular  frequency  to.  equals  the  angular  frequency 

li)   ..   Since  quite  accurate  common  frequency  signals  from  the 
K.+  X 

two  Omega  stations  are  necessary  for  navigation,  to,,.  .  .  =  to. 
will  be  assumed. 
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For  GRAN,  an  acquisition  and  reference  signal  (A/R  tone) 
will  be  generated  in  the  SARCOM  and  transmitted  with  the 
Omega  data.   The  transmitted  information  will  occupy  2«5  kHz, 
Figure  2,  a  smaller  bandwidth  than  the  3.4  kHz  received 
signal  bandwidth,  for  spectrum  conservation  and  will  be 
centered  at  a  frequency  in  the  allocated  band  of  frequencies 
from  406.0  -  406.1  MHz.   The  A/R  tone  will  serve  two  func- 
tions.  For  a  short  period  at  the  beginning  of  the  transmis- 
sion period  it  will  be  radiated  alone  by  the  SARCOM  so  the 
SARCEN  phase-locked-loop  receiver  can  achieve  lock  and  track 
the  subsequent  transmission.  The  other  function  of  the  A/R 
tone  is  to  allow  phase  shifts  and  frequency  drifts  in  the 
transmission  channel  through  the  SARSAT  to  the  SARCEN  to 
be  eliminated.   The  particular  contributors  to  these  effects 
will  be  satellite  link  delays  and  frequency  drift  in  the 
satellite  frequency  translator.   Proper  processing  at  the 
SARCEN  will  then  result  in  A.  values  in  (4)  being  caused 
only  by  relative  motion  of  the  SARCOM  with  respect  to  the 
Omega  stations.   The  frequency  drifts  in  the  SARCEN  or 
Omega  stations  during  the  transmission  period  will  be 
negligible  because  of  the  accurate  frequency  standards  which 
will  be  used. 

The  same  transmission  channel  will  be  used  to  receive, 
relay,  and  process  the  (K+i)th  frequency  as  the  ith  since 
they  must  be  at  the  same  frequency.   Therefore,  with  the 
exception  of  timing  errors  at  the  stations,  the  phase  shift 
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6     ■  9  ,  (i  =  1,  2 ,  K) .   The  phase  shifts  in  the  K 

different  channels  will  result  from  different  phase  shifts 
at  different  frequencies  in  the  transmission  channel  from 
the  SARCOM  to  the  SARCEN  and  from  the  frequency  division  in 
the  SARCOM.   The  compacting  of  the  spectrum  in  the  SARCOM 
must  be  accomplished  by  multiplication  in  each  receiving 
channel  with  an  appropriate  frequency  signal  created  by  fre- 
quency division  of  a  local  oscillator  signal.   The  frequency 
division  will  be  done  by  digital  counting  circuits  whose 
initial  counting  level  will  be  a  random  quantity.   The  phase 
shifts  in  the  K  channels  which  result  will  thus  be  mutually 
independent  random  variables,  uniformly  distributed  on  the 
interval  (-it,  it].   It  can  be  shown  that  this  assumption  is 
not  affected  by  the  addition  of  constant  phase  shifts  to  the 
signals  due  to  Omega  station  timing  errors  since  a  uniformly 
distributed  random  variable  plus  a  constant  remains  uniformly 
distributed . 

The  propagation  phase  shifts  of  the  signals  from  station 
one  and  station  two  will  be  written  as 


f 

'  .    =    to .      /  f—r- 

x  i    J  v±(y) 


-    2ttiu.  1    <    i    <    K       (5a) 


and 


* 


K+ 


1    J  VK+i(T1 


r-   -2Trm__.  .       1    <    i    <    K    .  (5b) 
)  K+l 


In  (5),  v.(  )  is  the  effective  phase  velocity  of  the  ith 
signal  from  the  Omega  station  to  the  SARCOM  receiver.   These 
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phase  velocities  are  written  as  a  function  of  the  distance 
from  the  transmitter  to  allow  for  incorporation  of  changes 
in  the  predicted  value  of  the  phase  velocities,  e.g.,  from 
daylight  path  to  dark  path  conditions.   The  integers  m.  are 
the  number  of  whole  cycles  of  phase  delay  the  ith  signal 
has  experienced  from  the  Omega  transmitters  to  the  SARCOM 
receiver.   The  values  x ..  and  x„  are  the  great  circle  dis- 
tances from  the  respective  Omega  stations  to  the  SARCOM. 

It  is  convenient  since  only  an  LOP  estimate  can  result 
from  two  stations,  to  estimate  only  one  variable  rather  than 
two  (x.  and  x„)  in  (5).   The  graphical  description  of  this 
method  is  contained  in  Chapter  IV.   It  consists  of  defining 
the  sum  x   +  x   =  L  and  replacing  x1  by  x.   Then,  by  holding 
the  value  of  L  constant  the  only  variable  to  estimate  is  x. 
Substitutions  of  these  definitions  into  (5)  yield 


x 

i     i  /     v±(Y) 


-2Trm.     1  _<  i  <    K   (6a) 


and 


^K+i  =  U 


^k 7  "27rnW-   1  <  i  <  K.  (6b) 

vK+i(L-y)      K+i     -    - 


Equations  (5)  and  (6)  are  written  as  two  equations  in  order 

to  emphasize  that  a) .  =  wT,  ,  .  and  to  allow  the  limits  of  the 
r  l     K+i 

integrals  to  be  written.   The  value  of  x  is  the  great  circle 
distance  from  the  first  station  to  the  receiver  at  the 
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beginning  of  the  data  collection  period  (t=0)  and  L  is  the 
sum  of  the  great  circle  distances  to  the  two  stations  for 
t  =  0 .   Any  SARCOM  motion  must  be  accounted  for  by  using  the 
doppler  estimation  information  from  Appendix  A  and  the  posi- 
tion est  ima  t e . 

The  large  range  over  which  VLF  signal  strengths  can  be 
expected  to  vary  and  the  impulsive  nature  of  the  atmospheric 
noise  at  VLF  frequencies  [11]  dictate  the  use  of  a  liraiter 
in  the  receiver.   This  thesis  will  assume  that  the  receiver 
will  use  a  hard  limiter  which  will  yield  the  necessary 
dynamic  range  as  well  as  minimize  the  effect  of  the  impulsive 
noise . 

The  components  of  the  noise  column  vector  N(t)  in  (3) 
will  be  assumed  to  be  sample  functions  from  white  Gaussian 
random  processes,  the  ith  component  of  N(t)  having  two-sided 
noise  spectral  density  height  of  N./2.   The  ith  and  (K+i)th 
noise  components  will  be  received  in  the  same  channel  of  the 
receiver    since  they  are  associated  with  signals  of  equal 
frequency.  But, as  shown  by  Figure  3, they  are  received  at 
different  times  and  will,  therefore  be  uncorrelat ed .   The 
noise  will  actually  be  the  sum  of  thermal  and  atmospheric 
noise  and  the  rationale  for  the  Gaussian  assumption  is  as 
follows.   The  atmospheric  noise,  expected  to  dominate  at  the 
VLF  frequencies,  will  be  the  sum  of  many  atmospherics  from 
around  the  globe.   The  distant  noise  impulses  will  have  in- 
sufficient amplitude  to  dominate  the  hard  limiter  in  the 
receiver,  but  nearby  impulses  will  dominate  the  hard  limiter. 
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Since  VLF  signals  propagate  very  well  over  long  distances, 
the  small  amplitude  impulses  of  noise  will  occur  quite  fre- 
quently and  by  the  central  limit  theorem  their  sum  will 
approach  a  Gaussian  random  process.   The  strong  impulses 
which  dominate  the  hard  limiter,  yielding  approximately  zero 
SNR  to  the  receiver  during  their  existence,  will  occur  much 
less  frequently  and  with  a  wise  choice  of  bandwidth  for  the 
K  channels  of  the  receiver,  each  impulse  will  capture  the 
hard  limiter  for  only  a  small  time.   It  can  be  shown  that 
this  time  will  be  on  the  order  of  the  inverse  of  the  band- 
width of  the  filter,  so  that  a  100  Hz  filter  will  allow  the 
hard  limiter  to  be  captured  for  a  time  of  about  0.01  sec. 
by  an  impulse. 

The  2K  dimensional  phase  velocity  vector  V  will  be  used 
to  represent  the  effective  phase  velocities  which  are 
modeled  as  jointly  Gaussian  random  variables.   The  covariance 

matrix  of  the  components  of  V  will  be  written  as  K  .   Since 

v 

inadequate  evidence  is  presently  available  to  justify  differ- 
ent uncertainties  for  the  effective  phase  velocities  when 
modeled  in  this  fashion,  it  is  convenient  for  the  simulations 
reported  on  in  Chapter  V  to  assume 


K   =  E  [V  V  ]  =  C 
v 


0 


0 


(7) 


In  (7),  a    is  the  scalar  variance,  assumed  equal  for  all 
signals.   The  symbol  p  represents  the  K-by-K  correlation 
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matrix  for  the  signals  from  a  single  station.   The  symbol  0 
represents  a  K-by-K  zero  matrix,  indicating  interstation 
correlation  is  negligible.   The  rationale  and  justification 
for  these  assumptions  is  contained  in  Chapter  I. 


42 


\ 


III.   AMBIGUITY  RESOLUTION 

Phase  differencing  systems  such  as  Omega  necessarily  are 
ambiguous.   The  size  of  the  unambiguous  lanes  (the  lanes 
formed  by  indistinguishable  LOP  from  a  station  pair)  is  deter- 
mined by  the  choice  of  frequencies.   A  discussion  of  how  the 
frequencies  determine  the  lane  size  is  contained  in  Chapter  I. 

The  Radux  navigation  system  {9]  which  immediately  pre- 
ceded the  Omega  system  had  a  200  Hz  modulation  of  its  carrier 
frequency,  40  kHz.   This  allowed  for  nominally  400  nmi  un- 
ambiguous lanes  which  were  considered  to  be  adequate  for  navi- 
gation.  It  was  the  short  baselines,  due  to  the  short 
propagation  ranges,  which  caused  Radux  to  be  abandoned  and 
not  its  ambiguities.   The  short  baselines  caused  the  diver- 
gence of  the  LOP  off  the  baseline  to  occur  very  rapidly  as 
distances  from  the  baseline  increased,  thus  reducing  the 
accuracy  of  the  system.   When  Omega  was  in  the  design  stage, 
dead  reckoning  (DR)  systems  and  the  reliability  of  receivers 
were  improving  markedly.   This  influenced  the  final  design 
of  Omega  to  concentrate  on  one  or  two  frequencies  which 
would  be  used  for  navigation.   Many  Omega  navigators  presently 
use  only  the  10.2  kHz  Omega  signal.   Also,  the  Defense  Mapping 
Agency  only  publishes  propagation  prediction  correction 
factors  for  the  10.2  kHz  and  13.6  kHz  Omega  signals.   As 
shown  in  Chapter  I,  the  three  navigation  signals  presently 
radiated  by  the  Omega  stations  yield  nominally  72  nmi  un- 
ambiguous lanes  at  the  baseline. 
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This  has  been  adequate  for  all  systems  using  Omega  to 
date.   The  SAR  applications  such  as  GRAN,  which  prompted 
this  study,  however,  require  an  ambiguity  resolution  tech- 
nique in  order  to  be  successful.   Since  the  retransmitter 
at  the  SAR  incident  will  in  general  have  no  prior  knowledge 
of  position  which  can  be  relayed  to  the  SAR  forces,  the 
assumption  must  be  made  that  complete  position  estimation 
must  be  accomplished  from  the  retransmitted  Omega  data. 
There  are  presently  five  methods  which  are  seriously  being 
investigated  as  potential  ways  to  solve  the  ambiguity 
problem.   They  are:   1)  addition  of  frequencies  to  the  Omega 
format,  2)  pulse  time-of -arrival  measurement,  3)  ratio  of 
signal-to-noise  ratios,  i.e.,  signal- to-s ignal  ratios,  4) 
multiple  LOP  examination,  and  5)  modulation  of  one  of  the 
Omega  signals  with  a  large  period  tone.   The  fifth  technique 
has  been  stated  by  the  Omega  Project  Office  in  Naval  Elec- 
tronic Systems  Command  as  being  unacceptable  to  the  host 
countries  for  the  Omega  stations  and  will  not  be  considered 
here  in  detail.   The  United  States  Air  Force,  however,  is 
funding  an  experiment  which  Cincinnati  Electronics  Corpora- 
tion and  NELC  of  San  Diego  are  to  carry  out  by  1975  which 
will  modulate  one  of  the  signals  from  the  experimental 
Omega  station  in  Forestport,  New  York  with  a  226  2/3  Hz  tone. 
Since  the  Omega  Project  Office  has  stated  that  a  final  imple- 
mentation using  modulation  on  an  Omega  signal  is  not  feasible, 
the  experiment  appears  to  be  academic. 
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As  discussed  in  Chapter  I ,  the  choice  of  frequencies 
determines  the  size  of  the  unambiguous  lanes.   It  is  reason- 
able then  to  expect  that  different  or  additional  frequency 
signals  for  Omega  could  yield  larger  unambiguous  lanes.   It 
is  very  likely,  however,  that  the  three  present  frequencies 
must  remain  and  that,  for  political  reasons,  at  most  only 
one  more  frequency  common  to  all  stations  can  be  added  to 
the  Omega  format.   This  statement  is  made  because  of  the 
many  existing  receivers  designed  for  the  three  existing 
Omega  frequencies  and  the  great  demand  for  the  remaining 
Omega  format  slots  for  other  applications.   Chapter  V  of 
this  dissertation  is  dedicated  to  the  study  of  the  perform- 
ance of  four-frequency  Omega  for  several  easily  implemented 
fourth  frequencies.   One  additional  frequency  is  not  ex- 
pected to  be  able  to  completely  solve  the  ambiguity  problem; 
however,  it  can  place  much  less  demand  on  other  ambiguity 
resolution  techniques  by  expanding  the  unambiguous  lane 
size  . 

Reference  [21]  discusses  the  use  of  the  relative  times- 
of-arrival  of  the  Omega  pulses  from  station  pairs  at  a  re- 
ceiver to  determine  a  lane  of  LOP  which  contains  the  receiver, 
One  of  the  conclusions  of  [21]  is  that  the  existing  three- 
frequency  Omega,  72  nmi,  baseline  ambiguity  cannot  be  re- 
solved in  high  noise  areas  with  .95  confidence  with  less  than 
on  the  order  of  one  hour  of  data.   The  pulse  time-of -arrival 
technique  is,  however,  expected  to  be  able  to  resolve  nominal 
360  nmi  baseline  lanes  with  three  minutes  of  Omega  data. 
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These  lanes  result  from  four-frequency  Omega  utilizing  10.88 
kHz  in  addition  to  the  existing  Omega  tones.   This  resolution 
is  expected  even  for  worst  case  noise. 

Reference  [22]  reports  on  an  experiment  which  examined 
the  potential  Omega  lane  size  which  could  be  resolved  using 
the  relative  powers  of  the  signals  received  from  a  station 
pair.   The  OPLE  ground  station  in  Dallas,  Texas  was  utilized 
in  measuring  satellite  relayed  Omega  data  retransmitted  by  a 
bread-boarded  SARCOM  from  various  locations  on  the  globe. 
It  was  assumed  that  the  noise  averaged  to  a  constant  over 
the  essentially  concurrent  measurement  interval  for  signals 
from  two  stations.   The  ground  station  was  used  to  estimate 
the  SNR  from  the  two  stations  and  with  the  above  assumption 
about  the  noise,  the  ratio  of  the  SNR's  from  the  stations 
for  a  given  frequency  is  a  signal  power  ratio.   The  conclu- 
sion of  [22]  was  that  signal- to-signal  power  ratios  are 
feasible  for  use  in  ambiguity  resolution,  but  the  limited 
experiment  performed  did  not  allow  for  a  conclusion  about 
the  lane  size  which  could  be  resolved.   The  SNR  for  the  ex- 
periment was  measured  at  the  ground  station  by  a  non-optimal 
technique  and  it  is  possible  that  the  maximum  likelihood 
estimate  of  SNR  derived  in  Appendix  B  of  this  dissertation 
could  be  used  to  advantage  if  this  technique  is  applied  in  a 
final  system. 

Reference  [23]  reports  on  simulations  which  were  carried 
out  for  three-frequency  multiple  station  Omega  utilizing 
multiple  LOP  intersections  to  resolve  ambiguity.   All  the 
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polygons  formed  by  LOP  from  stations  treated  pairwise  are 
examined  for  an  area  of  radius  250  nrai  on  the  globe  and  the 
one  which  maximizes  a  likelihood  function  is  selected  as 
being  the  most  likely  position  estimate.   The  likelihood 
function  is  formed  by  assuming  additive  Gaussian  phase  un- 
certainties affecting  the  phase  differences  which  correspond 
to  the  LOP  estimates  of  the  receiver  location  when  the  Omega 
stations  are  treated  pairwise.   The  conclusion  drawn  from 
the  simulations  is  that  the  multiple  LOP  method  of  ambiguity 
resolution  does  have  promise  in  terms  of  resolving  the  Omega 
ambiguities  if  multiple  position  estimates  along  with  their 
associated  relative  probabilities  are  taken  as  a  solution 
set . 

Additional  information  concerning  the  simulations  reported 
upon  in  [23]  is  contained  in  [24,  25].   It  can  be  found  from 
these  references  that  the  standard  deviations  used  for  the 
additive  Gaussian  phase  uncertainties  for  the  simulations 
was  0.003291  rad  for  three  stations  and  0.010408  rad  for  the 
remainder.   This  corresponds  to  LOP  estimate  standard  devia- 
tions on  the  order  of  0.0047  to  0.0142.   These  values  are 
said  to  be  derived  from  a  10  db-Hz  SNR  expected  from  three 
stations  and  0  db-Hz  SNR  from  the  remaining  stations.   It 
was  found  by  experiment  and  is  discussed  in  a  later  chapter 
of  this  thesis  that  standard  deviations  of  the  LOP  estimate 
on  the  order  of  0.1  to  0.2  baseline  nmi  are  experienced  with 
daylight  paths  and  SNR  in  excess  of  10  db-Hz  for  a  station 
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pair  using  the  optimum  estimator.   The  equation 


a.    =  0  (2tt/360) 
<p     x 


(8) 


allows  these  standard  deviations  to  be  converted  to  standard 
deviation  of  phase  difference  since  2 u  radians  of  phase  dif- 
ference corresponds  to  one  unambiguous  lane  on  the  baseline 
(nominally  360  nmi  for  the  four  frequencies  which  were  used 
in  the  experiment).   Therefore,  standard  deviations  of  phase 
difference  uncertainties  which  were  experienced  in  the  experi- 
ment were,  for  best  conditions,  on  the  order  of  0.002  to 
0.004  radians.   When  unambiguous  lane  size  is  increased  by 
the  addition  of  a  fourth  frequency  as  was  the  case  for  the 
experiment,  the  same  number  of  polygons  as  were  considered 
in  [23]  occur  in  a  larger  area  on  the  globe.   If  this  area 
is  expanded  sufficiently  both  by  addition  of  frequencies  and 
consideration  of  more  polygons,  the  multiple  LOP  method  of 
ambiguity  resolution  may  be  a  viable  technique  for  SAR  sys- 
tems . 

The  pulse  time-of -arrival  technique  and  the  multiple  LOP 
technique  appear  to  have  the  most  promise  as  ambiguity  reso- 
lution techniques  which  can  be  implemented.   The  multiple 
LOP  is  an  alternative,  although  nonoptimal,  technique  for 
treating  the  multiple  station  problem  to  be  discussed  in 
Chapter  VII.   The  conservation  of  computer  time  which  it 
allows  would  be  a  reason  for  adopting  it  as  a  method  over 
the  optimum  estimator.   An  estimation  procedure  using  the 
multiple  LOP  technique  would  be  to  obtain  the  two-station 
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LOP  estimates  optimally  and  determine  the  set  of  most  likely 
position  estimates  from  a  multiple  LOP  approach  as  described 
in  [23].   Use  of  pulse  time-of -arrival  information  with  this 
approach  will  be  very  helpful  in  reducing  the  area  on  the 
globe  over  which  LOP  intersections  must  be  examined. 
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IV.   TWO-STATION  LOP  ESTIMATION 

A.   TWO-STATION  ESTIMATION  FOR  K  FREQUENCIES 

In  order  to  establish  an  estimation  problem  which  will 
yield  an  unambiguous  estimate  of  LOP  for  two  stations  and 
K  frequencies,  we  assume  that  an  ambiguity  resolution  tech- 
nique has  resolved  the  lane  ambiguity.   Thus,  the  estimation 
is  performed  within  one  unambiguous  lane.   Within  the  pre- 
supposed unambiguous  lane  we  can  only  estimate  a  LOP  for  any 
station  pair.   Therefore,  referring  to  the  mathematical  mo- 
del of  Chapter  II,  an  essentially  equivalent  problem  is  to 
fix  L  in  (6a)and  vary  x  in  (6a)and  (6b)over  a  range  of  values 
necessary  to  traverse  all  candidate  LOP  in  the  unambiguous 
lane  given  for  the  two  stations  being  considered.   Figure  6 
illustrates  the  geometry  of  the  problem.   This  method  is 
convenient  for  the  following  reasons.   First,  only  one  para- 
meter, x,  need  be  varied  in  the  estimation  procedure;  and 
second,  for  a  fixed  step  size  in  an  iterative  solution,  the 
total  range  over  which  x  is  varied  (one  unambiguous  baseline 
lane  width)  is  the  same  for  any  receiver  position.   A  reason- 
able value  to  assign  to  L  is  the  sum  of  the  great  circle  dis- 
tances from  the  two  stations  to  the  center  of  the  presupposed 
area.   If  x   is  defined  as  the  great  circle  distance  from 
the  reference  station,  henceforth  referred  to  as  station  one, 
to  the  center  of  the  presupposed  area,  the  range  of  values 
over  which  x  must  be  varied  is  x   +  one-half  the  unambiguous 
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lane  size.   The  fact  that  this  range  of  values  for  x  actually 
traverses  an  unambiguous  lane  is  easily  seen  if  one  assumes 
two  artificial  stations  separated  by  a  baseline  distance  L. 
Estimation  of  x  on  the  artificial  baseline  formed  is  equiva- 
lent to  estimating  x  off  the  baseline  in  the  real  problem  if 
the  manner  of  estimation  explained  above  is  followed.   Thus, 
the  range  of  x  which  must  be  examined  is  one  baseline  unam- 
biguous lane  width. 

To  estimate  x  in  the  problem  posed,  we  form  the  likeli- 
ho6d  function  [26]  as 


A(x)  =  f(Y/x) 


// 


f (Y/V,0,x)f (V)f (G)dVd6 


(9) 


where  V  is  the  2K-dimens ional  phase  velocity  column  vector, 
0  is  the  K  dimensional  random  phase  vector  and  f(V)  and  f(6) 
are  the  joint  probability  density  functions  of  the  components 
of  V  and  6  which  are  assumed  mutually  independent.   f (Y/x)  is 
the  joint  probability  density  of  the  components  of  the  vector 
Y  conditioned  upon  knowing  x.   Since  our  only  interest  in 
the  likelihood  function  is  to  find  the  values  of  x  which 
yield  its  maxima,  we  can  equivalently  maximize  any  monotonic 
function  of  it.   Factors  multiplying  the  likelihood  function 
but  not  containing  x,  v.  or  6 .  will  be  dropped  without  re- 
defining the  likelihood  function  throughout  the  dissertation. 

Ambiguity  resolution  information  can  be  conveniently  in- 
corporated in  the  maximum  likelihood  estimator  in  the  follow- 
ing way.   Consider  the  ambiguity  resolution  information  as 
a  priori  information  about  x  in  the  form  of  a  probability 
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density  f(x).   Equation  (9)  when  multiplied  by  f(x)  is  pro- 
portional to  the  a  posteriori  density  of  x.   The  maximum  a 
posteriori  (MAP)  estimate  of  x  is  then  that  x  which  maxi- 
mizes this  equation  and  is  the  optimum  estimate  of  x  for 
the  problem  [26].   When  f(x)  is  a  uniform  density  function, 
as  would  be  the  case  if  exactly  one  unambiguous  lane  were 
given  by  an  ambiguity  resolution  technique  with  no  informa- 
tion about  where  in  the  lane  the  correct  LOP  might  be,  the 
MAP  estimator  is  equivalent  to  the  maximum  likelihood  esti- 
mator . 

It  can  be  shown  [26]  that  for  the  density  function  of 
N(t)  assumed  for  the  mathematical  model  in  Chapter  II, 

T 
2K  f 

f(Y/V,6,x)oc  1F   exp[(2/Ni)  /  y±(t) s±(t ,x)dt]        (10) 

i=l  J 

o 

where  T  is  the  total  data  collection  period.   Substituting 
(4)  in  (10)  and  using  the  trigonometric  identity  for  the 
cosine  of  the  sum  of  angles,  (10)  can  be  written  as 

2K       , 
f(Y/V,6,x)cc  n   exp  <(2Ai/Ni)[Ci  cos  ((JK-S^ 
i=l 


+  S.  sin  ((f). -0  .)]  \ 
x       ri   1   ) 


(ID 


In  (11)  C.  and  S.  are  defined  as 


M 
j  =  l 


10(j-l)  +  T. 

y.(t)  cos[(u.+A.)t]dt    (12a) 
10(j-l) 
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and 


10(j-l)  +  Tx 

y  At)    sin  [  (w.+A  .  )  t  ]dt 
10(j-l) 


(12b) 


where  M  is  the  number  of  T.  sec.  duration  pulses  contained 

1 

in  the  data  collection  period. 

The  quantities  A.,  N.  and  A.  (l<i<2K)  must  be  estimated 

from  the  collected  data.   The  estimation  equations  for  these 

quantities  are  derived  in  Appendix  A  and  Appendix  B.   From 

this  point  in  the  two-station  estimation  derivation  it  will 

be  assumed  that  the  values  of  A.  and  N.  have  been  estimated 

l       i 

from  the  data  and  the  estimates  are  available.   The  values 

of  A.  will  be  assumed  to  have  been  estimated  and  the  data 
l 

corrected  for  doppler  shift  utilizing  this  information.   Fol- 
lowing these  assumptions,  we  define 


R±2  =  (C.2  +  S.2)(2A./N.)2 


Ki<2K 


(13a) 


and 


u.  =  tan_1(S./C.) 


Ki<2K. 


(13b) 


The  likelihood  function  (9)  may  then  be  written  by  substitu- 
tion of  (13)  into  (11)  into  (9)  as 

K 


A(x)    =    Ev/f(8)exp|    J^     [Ricos(<J)i-ui)+Ri+Kcos(cf>i+K-yi+K)]cosei 

i=l 

K  , 

+  V*     [R.sin(((i.-y.)+R..1,sin(((i..„-v..ir)]sin    0.>d9 
Z-/  i  i       i  i+K  ri+K       l+K  l) 

i=l 

(14) 


where  E  (•)  denotes  expectation  with  respect  to  V. 
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Since  the  9   were  argued  to  be  independent,  uniformly 
distributed  random  variables  on  the  interval  (-7T,  Tr]  in 
Chapter  II,  we  may  use  the  relationship 

IT 

(2tt)-1   /  exp[A  cos  9  +  B  sin  9]d9  =  IQ  [  (A2+B2 ) l!  2  ]  ,  (15) 

-ir 
some  algebra  and  some  trigonometric  identities  to  write  the 

likelihood  function  as 

.     A(x)  =  Ev   5   Io{[R2  +  R2+i+  2R.RK+icos  *.]1/2}      (16) 
i=l 


where 

\l>  .    =  principal  value   of  [^K+i_VJ;L-(i>K+:L+*i  1  •         (17) 

In  (15)  and  (16)  I  (•)  is  the  zero  order  modified  Bessel 
function  of  the  first  kind.   At  this  point,  no  approximations 
have  been  made  and  the  assumed  random  nature  of  the  effec- 
tive phase  velocities  has  not  been  utilized.   Therefore,  if 
the  phase  velocity  column  vector  V  is  assumed  to  have  deter- 
ministic components  as  was  first  thought  to  adequately  model 
the  problem,  the  maximum  likelihood  estimation  equation  for 
known  phase  velocities  can  be  obtained.   This  estimation 
equation  is  derived  in  Appendix  C.   As  was  discussed  in 
Chapter  I,  the  phase  velocities  are  not  well  enough  known 


1The  value  of  9  which  lies  in  the  interval  -tt<9_<tt  is 
called  the  principal  value  of  9. 
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to  use  this  model,  so  we  continue  with  the  derivation  of  the 
estimator  assuming  Gaussian  phase  velocities. 

A  few  minutes  of  data  collection  is  expected  to  be  ade- 
quate to  insure  high  SNR  in  one  Hz  even  for  worst  case  noise 
[2],  The  arguments  of  the  I  function  associated  with  the 
ith  frequency  are  on  the  order  of  the  sum  of  the  SNR  in  one 
Hz  of  the  signals  received  from  the  two  stations  at  the  ith 
frequency.  This  allows  the  replacement  of  the  I  functions 
in  (16)  by  their  large  argument  approximation  [27], 


I  (z)  -  (2ttz)  1/2  exp(z). 
o 


(18) 


Using  (18),  the  likelihood  function  (16)  may  then  be  written 


as 


exP{  E  [R2  +  R2+.  +  2R.RK+.cos  ^.]1/2} 

A(x)  =  Ev   -   1  =  1 -(19) 

V   K    ?  2  1/4 

J^i  +  RK+i  +  2RiRK+i  C°S  *1] 


The  numerator  of  (19)  is  a  rapidly  varying  function  of  x, 
sharply  peaked  at  all  the  local  maxima  of  the  function.   The 
denominator  is  a  more  slowly  varying  function  and  for  regions 
near  a  local  maximum  can  be  considered  to  be  nearly  constant. 
The  denominator  will  have  different  values  from  peak  to  peak 
of  the  likelihood  function,  but  the  magnitude  of  the  change 
from  peak  to  peak  will  be  quite  small  compared  to  the  change 
in  the  magnitude  of  the  numerator  for  the  same  peaks.   For 
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the  peaks  in  the  vicinity  of  the  global  maximum  which  we  are 

seeking  to  find,  the  denominator  in  (19)  can  be  considered 

to  have  relatively  little  effect  on  the  likelihood  function 
and  will  be  ignored,  yielding 


K 


1/2 


A(x)  =  Ev  exp{  I     [R:    +  RK+±  +  2R±RK+i  cos  ^    >•   (20) 


i  =  l 


It  is  convenient  to  write 


v±(x)  =  v±(x)  +  pi 


(21) 


where  v  .  (x)  is  the  expected  value  and  p.  some  random  pertur- 
bation of  the  ith  phase  velocity.   The  expected  value  of  the 
phase  velocities  are  all  very  near  c  [10],  so  we  can  approxi- 
mate (6a)  and  (6b)  by 


x 


dY 


i(y)[l+£i] 


-  27rm. 

l 


Ki<K 


L 
.     ~   u      [   ^ 

K+i "  V  wl-^)[i+w 


(22a) 


-  2ffmv..    Ki<K     (22b) 


X 


where 


ei  =  pi/c  S  Pi^U) 


(23) 


The  random  variables  £.,  l£i£2K  will  be  taken  as  jointly 
Gaussian  with  zero  means  and  as  components  of  the  2K  dimen- 
sional vector  e.   This  follows  directly  from  the  assumption 
of  the  components  of  V  being  jointly  Gaussian  and  the 
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definition  of  e.   The  covariance  matrix  of  e  will  be 


K 


T  Kv 


(24) 


where  K   is  the  covariance  matrix  of  the  2K  dimensional 
v 

vector  V.   Utilizing  this  result  and  the  approximation 

— —  **  1  -  n  for  |n|<<l,  we  can  now  write 
l+n  '  '    ' 


0)  .x 
Yi    ri     c    l       i 


Ki<K 


(25a) 


and 


0).  (L-x) 

*K+i  "  ?K+i c £K+i  -  2™K+i   ^i* 


(25b) 


where 


_dx_ 
v±(y) 


Ki<K 


(26a) 


and 


K+ 


.  =  OK  fL *I 


Ki<K, 


(26b) 


Recalling  the  definition  of  if),  from  (17),  for  high  SNR 

if;  .  will  be  small  near  the  local  maxima  of  the  likelihood 
l 

function  and  we  write 


COS  if)  .  -  1  -  —j-  . 


(27) 
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This  allows  the  likelihood  function  to  be  written  as 


A(x)  =  E£  exp  {  E  (R±  +  RR+i)[l- 


1  K+1  \]1/2}.     (28) 
(Ri+RK+i) 


Again,  using  the  fact  that  ty  .    and  thus  \p  .     is  small  in  the 
vicinity  of  the  local  maxima  of  the  likelihood  function  and 
the  approximation 


1/2  ~ 

(l  -  n)  '     -  i  -  n/2, 


(29) 


an  .  equivalent  likelihood  function  can  be  written  as 

K       2 
A(x)  =  E   exp  [  E   Q,<K  ] 

fc       i=l 


(30) 


where 


Q±  -  - 


R  .  R,_ .  . 

l  K+i 


2(R.+R^.) 
l   K+ 1 


(31) 


Using  vector  notation,  the  likelihood  function  can  be 
written  as 


A(x)  =  E   exp  [e  Ae  +  BTe  +  D]. 


(32) 


The  2K-by-2K  matrix  A  is  defined  by  j th  row  and  nth  column 

entries  _   _ 

to  .   x  Q  . 

-2— !  1  <  j  =  n  <  K       (33a) 


jn 


'jn 


2     2n 

_  yK*  Q.j-K 

2 
c 


K+1  £  j  =  n  <    2K      (33b) 
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a 


2 

-0).   x(L-x)Q. 
-1  3 


Jn  2 


1  <  j  =| (n+K) |  <  2K  (33c) 


a,   =   0  elsewhere.        (33d) 


The  2K-dimensional  column  -vector  B  is  defined  by  ith  entry 


b,    =    -   -   Q.W.co.x  1    <    i    <    K  (34a) 

icxiii  — — 


bi    =    c    Qi-KWi-KWi-K(L_x)       K+1   -   i   -    2K  (34b) 


where 


W.  =  principal  value  of  [y  .  .  -y  . -d>  .  .  T.+(b  .  ]  .  (35) 

x         t-  t-  L  "l+K   x  Tx+K  ri 


The  scalar  D  is  defined  as 


K 
D  =   Z   Q.  W.  .  (36) 

i=l 


The  expected  value  in  (32)  under  the  conditions  given  is 
shown  in  Appendix  D  to  yield 

exp[D+(l/2)BT  K  (I-  2AK  )-1B] 
A(x)  =  £ £ (37) 

|I  -  2AK£|1/2 

where  I  is  the  2K  by  2K  identity  matrix  and  |*|  denotes  the 
determinant  of  a  matrix.   The  value  of  x  which  maximizes 
the  likelihood  equation  can  be  determined  by  allowing  x  in 
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(37)  to  range  over  the  necessary  values  as  shown  in  Figure  6. 
Other  less  likely  LOP  can  be  determined  and  a  measure  of  the 
probability  of  their  being  correct  can  be  found  by  dividing 
the  value  of  A  for  the  less  likely  LOP  by  the  value  of  A  for 
the  most  likely  LOP.   If  the  LOP  being  estimated  is  considered 
to  be  a  random  variable,  uniformly  distributed  across  one  un- 
ambiguous lane,  this  ratio  is  a  ratio  of  the  a  posteriori 
.probability  densities  of  the  two  LOP. 

Figure  7  is  a  drawing  of  the  envelope  of  a  typical  likeli- 
hood function  as  it  varies  with  x.   Within  the  "envelope" 
the  function  will  be  a  more  rapidly  varying  function  of  x. 
The  number  of  local  maxima  in  an  unambiguous  lane,  determined 
by  the  values  and  number  of  Omega  frequencies,  will  in 
general  be  many  more  than  shown.   The  value  of  x  ,  the  dis- 
tance from  the  first  station  to  the  center  of  some  presup- 
posed area,  serves  only  as  a  means  to  establish  the  fixed 
value  of  L  to  use  in  the  estimation.   In  Figure  7,  x  .,  the 
value  of  x  which  yields  the  maximum  value  of  A  cannot  be 
used  as  a  valid  measure  of  the  distance  of  the  receiver  from 
the  reference  station,  but  serves  only  to  identify  the  re- 
sulting LOP  estimate. 

B.   TWO-STATION  ESTIMATION  FOR  FOUR  FREQUENCIES 

The  general  estimation  equation  for  K  frequencies  re- 
ceived from  two  stations  is  given  by  (37).   This  equation 
can  be  applied  to  the  four-frequency  case  with  the  size  of 
the  matrices  A  and  K   being  8-by-8  and  the  column  vector  B 
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being  8  dimensional.   For  the  four-frequency,  two-station 
problem  then,  the  likelihood  function  is  given  by 


A(x)  = 


exp[D+(l/2)BTK  (I-2AK  )  ^~    B] 


I  I  -  2AK  | 


1/2 


(38) 


where  K  is  the  8-by-8  covariance  matrix  of  the  normalized 
effective  phase  velocity  vector  £.  The  8-by-8  matrix  A  is 
given  by 
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where  Q  ,  1  <  ±  <_    4,  is  given  by  (31).   The  column  vector  B 
in  (38)  is  given  by 

~QlWla3lX 
-Q2W2w2x 

-^3W3W3X 

-Q4w4o)4x 


B  -* 

c 


Q1W10)1(L-x) 

^2W2W2  ^L_x^ 
Q3W3to3(L-x) 

Q4W4<VL-X) 


(40) 


where  W  ,  1  <  i  <  4,  is  given  by  (35).   The  scalar  D  is 
given  by 


D 


£   Q.W. 
i=l 


(41) 


C.   ESTIMATION  ERROR  TYPES 

The  cyclic  nature  of  the  likelihood  function  yields  LOP 
estimates  which  fall  into  three  categories: 

1.  Estimates  which  fall  within  nominally  1  or  2  nmi  or 
less  of  the  correct  LOP, 

2.  Estimates  which  are  nominally  one  half  wavelength 

of  the  fundamental  four  frequencies  from  the  correct 
LOP,  referred  to  as  "minor  lane"  errors  and 

3.  Estimates  which  are  farther  from  the  correct  LOP  than 
the  previous  two  types,  which  we  refer  to  as  "major 
lane"  errors. 
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The  latter  two  types  of  estimates  are  considered  as  erroneous 
estimates  for  the  purposes  of  this  thesis.   An  important  con- 
sideration to  note  is  that  the  probability  of  making  a  posi- 
tion error  when  using  data  from  several  stations  will  very 
likely  be  less  than  the  probability  of  making  a  lane  error 
in  the  two  station  case  for  the  following  reasons.   The  trans- 
missions of  at  least  three  Omega  stations  will  be  available 
(this  number  being  required  for  a  position  fix)  and  it  is 
very  likely  that  one  or  more  additional  stations  will  be 
usable  given  the  great  distances  of  VLF  propagation  and  the 
planned  worldwide  location  of  eight  transmitters  in  the  Omega 
system.   This  redundancy  will  reduce  the  probability  of  large 
position  estimation  errors.   The  multiple  station  problem  is 
considered  in  Chapter  VII. 
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V.   FOUR-FREQUENCY  SIMULATION 

In  the  interest  of  achieving  unambiguous  lane  widths 
which  can  be  identified  by  an  ambiguity  resolution  technique, 
additional  frequencies  may  be  added  to  the  Omega  format. 
Assuming  that  the  present  three  Omega  frequencies  are  re- 
tained, adequate  widening  for  ambiguity  resolution  by  time- 

f-arrival  techniques  [21]  is  expected  to  be  achieved  by 
addition  of  a  fourth  frequency.   A  fifth  frequency  and/or 
frequency  modulation  of  one  of  the  signals  by  a  signal  with 
wavelength  on  the  order  of  several  thousand  miles  would 
yield  better  estimator  performance  and  perhaps  inherent 
ambiguity  resolution  as  well.   This  chapter  considers  the 
performance  of  the  four-frequency  Omega  maximum  likelihood 
estimator  for  several  candidate  fourth  frequencies  which 
can  be  added  to  the  Omega  format  with  minimum  modification. 

A.   FREQUENCY  CHOICE 

The  criteria  used  in  selection  of  the  fourth  frequencies 

were: 

1.  Retention  of  the  present  three  Omega  frequencies 

(10.2  kHz,  11  1/3  kHz,  and  13.6  kHz), 

2.  Separation  of  the  fourth  frequency  from  any  of  the 
existing  three  Omega  frequencies  or  other  known 
existing  or  planned  VLF  transmissions  by  at  least 
250  kHz  to  avoid  interference  in  present  receivers, 
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3.  A  resultant  unambiguous  lane  size  greater  than  the 
present  three  frequency  baseline  size  of  nominally 
72  nmi  and 

4.  Selection  of  a  fourth  frequency  which  can  be 
generated  using  the  present  Omega  signal  format 
generator  frequency  of  816  kHz. 

In  addition  to  these  criteria,  preliminary  four-frequency 
simulations  were  carried  out  for  SNR  which  are  to  be 
expected  for  reasonably  long  transmission  paths  and  assum- 
ing that  exact  phase  velocity  predictions  were  available. 
These  simulations,  discussed  in  Appendix  C,  indicate  that 
four-frequency  unambiguous  lane  widths  must  be  restricted 
to  about  1000  nmi  or  less  for  estimation  performance  "accept- 
able" for  SAR  application.   The  frequencies  which  met  the 
four  criteria  above,  the  integer  divisor  of  the  fundamental 
station  frequency  (816  kHz)  and  the  resultant  nominal  base- 
line unambiguous  lane  width  are  presented  in  Table  III. 
Figure  8  shows  the  Omega  frequency  spectrum  with  the  present 
three  frequencies  and  the  potential  fourth  frequencies. 

B.   SIMULATION  DATA  AND  PROCESSOR 

Monte  Carlo  simulation  of  the  four-frequency  Omega  prob- 
lem was  carried  out  utilizing  the  maximum  likelihood  esti- 
mator given  by  (37)  in  Chapter  IV  with  K  =  4 .   The  parameters 
which  were  varied  were  the  choice  of  fourth  frequency,  the 
SNR  in  one  Hz  and  the  uncertainties  of  the  phase  velocities. 

It  is  shown  in  Appendix  A  that  C.  and  S.,  when  conditioned 

i       x 
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TABLE  III 
CANDIDATE  FOURTH  FREQUENCIES 

Nominal  Unambig 

Frequency  (f  <)              Divisor  (n.)  Lane  (1/2  A  ) 

10.4615  kHz                    78  928  nmi 

10.8800  kHz                    75  357  nmi 

11.6571  kHz                    70  500  nmi 
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upon  knowing  <j)   and  8   and  also  under  the  assumption  that 

N(t)  is  a  vector  whose  components  are  sample  functions  of 

a  white,  Gaussian  process,  are  Gaussian  random  variables. 

The  conditional  means  and  variances  of  C.  and  S.  when 

i       1 

doppler  shift  is  assumed  to  be  negligible  are  given  by 


A  MT. 
E[Ci|<J)i,ei]  =    2  X  cos  (-<}>.  +  8±)  ,  (42a) 

A.MT. 
E[S±|(J)i,ei]  =  -iy-i  sin  (-({k  +  0±)  ,  (42b) 

arid 

N.MT. 
VAR  [C.|c))1,9i]  =  VAR  [S1|$±,  6±]  =   X 4  X  .  (42c) 


For  the  simulation,  A.  and  9.  were  assumed  equal  to  zero  and 

A.  and  N.  were  determined  from  the  SNR  in  one  Hz,  assumed 
X        x  ' 

the  same  for  all  i.   The  product  MT .  was  assumed  equal  to 
18  for  all  i,  this  choice  being  influenced  by  the  GRAN  anti- 
cipated 180  second  data  retransmission  time,  equivalent  to 
18  nominally  one  second  pulses  of  Omega  data  for  each  fre- 
quency from  each  station.  With  these  assumptions,  the  means 

and  variances  of  C.  and  S.  conditioned  upon  <j>  .  and  6.  become 

xx  r    rx       x 

E  [C.U.,  6.]  =  9A.  cos  <J>.  ,  (43a) 

x'Txx       x      Tx' 

E  [S.U.,  6.]  =  9A..  sin  cf>  .  ,  (43b) 

x,Txx       1      rx 

and 

9N. 

VAR  [C.U.,  0.]  =  VAR[S.U.,  6.]  =  -^  .  (43c) 
x'x'   x  x|Tx    x      2 

The  eight  phase  velocities  used  to  generate  the  data  were 
generated  by  a  correlated  Gaussian  random  number  generating 
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computer  program  using  the  mean  vector  with  all  components 
equal  to  c  and  the  covariance  matrix 


2  2 

K  =  c  a 
v 


(44) 


In  (44)  p  is  a  symmetric  4  by  4  correlation  matrix  referring 
to  the  inter frequency  correlation  of  transmissions  from  each 

station  and  o  represents  a  4  by  4  zero  matrix  indicating  un- 

2 
correlated  interstation  phase  velocities.   In  (44)  a   is  the 

variance  of  all  e .  (equal  to  the  variance  of  v./c)  which  is 
varied  as  a  parameter  in  the  simulation.   The  correlation 
matrix  was  chosen  to  reflect  uncorrelated  interstation  trans- 
missions because  in  the  general  case  very  low  correlation  is 
expected  [19],   The  inter-frequency  correlation  matrix  p  was 
formed  by  assuming  the  correlation  of  10.2  kHz  and  13.6  kHz 
frequencies  from  one  station  to  be  .7,  which  seems  to  be  a 
reasonably  good  choice  [16,  18,  and  19]  and  linear  interpola- 
tion using  the  equation 


ij 


1  - 


.3 


3400 


(«±  -  v 


(45) 


where  p. .  is  the  ith  row,  jth  column  entry  in  the  matrix  p 
and  f.  is  the  ith  frequency.   These  assumptions  resulted  in 
the  symmetric  matrix 
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1.00 

.90 

.70 

.90 

1.00 

.80 

.70 

.80 

1.00 

P/.i 

P/.0 

Px* 

'14 


'24 


M34 
1.00 


(46) 


where  the  first  three  rows/columns  of  the  matrix  correspond 
to  the  frequencies  10.2  kHz,  11  1/3  kHz,  and  13.6  kHz  res- 
pectively.  The  values  used  for  p.,  =  P,.,  i  =  1,  2,  3  were 
chosen  from  Table  IV  depending  upon  the  fourth  frequency  be- 
ing used. 

The  values  for  L  and  x  used  in  the  creation  of  the  simu- 
lation data  were  8000  nmi  and  3000  nmi  respectively.   In  the 
final  Omega  system  implementation,  these  values  should  be 
typical  of  values  to  be  expected.   With  all  the  above  assump- 
tions, the  simulation  data  was  generated  by  computer  using 
random  number  generating  programs  for  the  random  variable 
sample  values. 

The  processing  of  the  data  with  the  four-frequency  maxi- 
mum likelihood  estimator,  was  accomplished  using  c  as  the 
mean  value  for  all  phase  velocities  and  the  covariance  matrix 


K 


T  Kv 


(47) 


where  K   is  given  by  (44)  and  the  value  of  the  parameter  a 
Using  c  for  the  mean  value  of  phase  velocity  for  processing 
is  equivalent  to  having  a  reasonably  good  predicted  mean 
value  of  phase  velocity  since  the  data  were  created  using 
this  mean  value. 

73  (   \ 


TABLE  IV 
FOURTH  FREQUENCY  CORRELATION  COEFFICIENTS 

Frequency  (f ,)  i  p 


14 


10.4615  kHz                    1  .98 

2  .92 

3  .72 

10.8800  kHz                    1  .94 

2  .96 

3  .76 

11.6571  kHz                    1  .87 

2  .97 

3  .83 
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Figures  9,  10,  and  11  illustrate  the  results  of  the  simu- 
lation for  each  of  the  three  choices  of  fourth  frequency. 
Each  point  plotted  consists  of  a  number  of  "runs"  which  cor- 
respond to  a  three  minute  data  collection  period.   The  frac- 
tion of  the  total  number  of  runs  for  each  point  which  yielded 
lane  errors  (either  major  or  minor,  as  defined  in  Chapter  IV) 
are  plotted  as  a  function  of  SNR  in  one  Hz  and  a,  the  standard 
deviation  of  the  phase  velocities.   Table  V  gives  the  number 
of  runs  which  were  utilized  in  the  generation  of  each  point 
and  a  breakdown  of  the  number  of  the  lane  errors  which  were 
major  and  minor. 

The  number  of  runs  which  were  used  in  generating  the 
points  in  Figures  9,  10,  and  11  are  quite  small  for  simula- 
tion, particularly  considering  the  small  number  of  resultant 
lane  errors.   The  reason  for  this  is  the  computer  time  which 
is  required,  15  to  40  seconds  on  the  IBM  360  per  run  depend- 
ing on  the  lane  size.   However,  the  reasonably  consistent 
increase  in  lane  errors  with  increasing  lane  width  and  phase 
velocity  standard  deviation  influences  the  decision  that 
the  general  results  are  valid. 

The  results  of  this  simulation  prompt  the  conclusion 
that  with  adequate  knowledge  of  the  phase  velocity  variance 
and  correlation  and  with  reasonable  signal-to-noise  ratios, 
the  unambiguous  lane  size  can  be  increased  beyond  the  present 
nominal  72  nmi  obtained  with  three  frequency  Omega.   The 
controlling  factors  in  the  selection  of  a  fourth  frequency 
for  Omega  to  be  used  for  SAR  purposes  are  (1)  the  capability 
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TABLE  V 

SIMULATION  RESULTS 

Lane  Errors 
Frequency  (f/)   Std.  Dev.(a)   SNR     No.  Runs   (ma j or ) (minor ) 

10.4615  kHz 


10.8800  kHz 


11.6571  kHz 
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of  ambiguity  resolution  techniques  to  resolve  the  lane  am- 
biguity with  restricted  amounts  of  Omega  data  and  (2)  the 
accuracy  with  which  the  phase  velocities  of  the  Omega  fre- 
quencies can  be  predicted. 

One  of  the  results  of  the  Monte  Carlo  simulation  as 
shown  by  Figures  9,  10,  and  11  is  that  as  unambiguous  lane 
size  is  increased,  the  probability  of  lane  error  increases. 
However,  the  estimator  provides  several  estimates  along  with 
their  relative  probabilities  of  being  correct.   Thus,  even 
when  a  lane  error  occurs,  the  second  or  third  most  likely 
estimate  is  nearly  always  correct.   This  fact,  in  conjunc- 
tion with  the  reduced  probability  of  position  error  expected 
when  redundant  information  from  additional  stations  is  avail- 
able, makes  a  .1  or  .2  probability  of  lane  error  more 
acceptable  as  a  trade-off  for  the  larger  unambiguous  lane 
size . 
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VI.   A  FOUR-FREQUENCY  OMEGA  EXPERIMENT 

All  of  the  ambiguity  resolution  techniques  discussed  in 
Chapter  III,  which  appear  practical  to  implement,  require  an 
unambiguous  lane  width  in  excess  of  the  nominal  72  nmi  which 
the  present  three  Omega  frequencies  yield.   The  performance 
which  can  be  expected  from  a  maximum  likelihood  estimator  if 
used  with  four-frequency  Omega  for  several  easily  implemented 
fourth  frequencies  was  presented  in  Chapter  V.   Naval  Air 
Test  Center,  Patuxent  River,  Maryland  conducted  a  four-fre- 
quency experiment  utilizing  one  of  these  frequencies  (10.88 
kHz)  in  the  fall  of  1973.   The  maximum  likelihood  estimator 
derived  in  Chapter  IV  was  utilized  to  process  this  data,  the 
results  of  which  are  reported  in  this  chapter. 

A.   DATA  COLLECTION 

The  four-frequency  Omega  data  was  collected  by  Texas  In- 
struments Corporation  at  Dallas,  Texas  [29].   The  receiver 
location  was  1197.53  nmi  from  the  Forestport,  New  York  experi- 
mental Omega  station  and  2349.62  nmi  from  the  Trinidad  Omega 
station.   The  data  was  collected  over  the  two  five-day  periods 
17-21  October  1973  and  14-18  November  1973.   Data  was  taken 
for  three-minute  intervals  four  times  per  hour  beginning  on 
the  hour  and  half  hour  and  seven  minutes  past  the  hour  and 
half  hour  for  the  entire  periods.   Each  of  these  three-minute 
collection  periods,  herein  referred  to  as  "runs,"  consisted 
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of  18  nominally  one-second  pulses  of  signal  for  each  frequen- 
cy from  each  station.   On  the  last  day  of  the  first  data  col- 
lection period,  the  Forestport  station  reduced  transmitted 
power  to  very  low  levels.   This  portion  of  the  data  is  not 
included  in  the  estimation  results  as  it  is  considered  not  to 
be  representative  of  actual  conditions  to  be  expected. 

The  receiver  used  for  the  data  collection  was  the  Omega 
Position  Locating  Equipment  (OPLE)  ground  station  modified  to 
receive  four  frequencies  [29].   A  block  diagram  of  the  receiv- 
ing equipment  is  shown  in  Figure  12. 

B.   MAXIMUM  LIKELIHOOD  ESTIMATION 

The  maximum  likelihood  estimator  for  high  SNR  Omega  was 
derived  in  Chapter  IV.   Quite  high  SNR  were  experienced  as 
shown  in  Appendix  B.   As  described  in  Chapter  II,  the  effec- 
tive phase  velocities  are  modeled  as  being  jointly  Gaussian. 
The  eight  element  column  vector  V  is  defined  as  the  effective 
phase  velocity  vector 

VT  =  [vx,  v2,  .  .  .  vg]  .  (48) 

The  first  four  elements  of  V  correspond  to  the  four  signals 
at  frequencies  10.2  kHz,  11  1/3  kHz,  13.6  kHz,  and  10.88  kHz, 
in  that  order,  from  Forestport.   The  last  four  elements  of  V 
correspond  to  the  same  frequency  signals  from  Trinidad.   The 
mean  values  of  the  Gaussian  vector  V  which  were  used  in  the 
estimation  were  those  provided  by  Professor  J.  A.  Pierce  of 
Harvard  in  [16]  and  in  private  correspondence.   These  values 
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are  given  as  ratios  of  the  phase  velocity  and  the  velocity 
of  light  in  vacuum  (v/c)  in  Table  II.   Because  of  the  signi- 
ficant diurnal  effects  of  phase  velocity,  different  values 
for  v/c  are  required  for  day  (all  lighted  transmission  path) 
and  night  (all  dark  transmission  path).   For  transition  per- 
iods, when  a  path  was  not  all  lighted  or  all  dark,  weighted 
averages  of  the  day  and  night  path  phase  velocity  values  from 
Table  II  based  upon  the  fraction  of  the  path  which  was  light 
was  used.   In  general,  this  would  not  be  justified,  but  since 
the  values  of  phase  velocity  for  day  and  night  are  so  close 
and  because  of  the  phase  velocity  uncertainties,  it  is  con- 
sidered acceptable  to  use  this  technique. 

Since  the  phase  velocities  are  considered  to  be  jointly 
Gaussian,  a  covariance  matrix  must  be  specified  in  order  to 
do  the  estimation.   The  inter-station  correlation  of  phase 
velocities  was  considered  to  be  zero  and  the  inter-frequency 
correlation  was  determined  by  assuming  the  phase  velocities 
of  the  13.6  kHz  signal  and  the  10.2  kHz  signal  to  have  a  0.7 
correlation  coefficient.   These  assumptions,  taken  from  [16, 
18,  and  19],  allow  the  matrix  p  to  be  determined  as  in  Chap- 
ter V,  Equation  (46) 
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The  matrix  p  Is  used,  as  in  Chapter  V,  to  form 


K 


2   2 
c   a 


(50) 


where  c  is  the  speed  of  light  in  vacuum,  O    is  the  standard 
deviation  of  v/c  which  was  assumed  to  be  the  same  for  all 
frequencies,  and  o  is  the  4-by-4  zero  matrix.   The  final 
parameter  which  must  be  specified  for  the  estimation  pro- 
cedure is  a.   As  discussed  throughout  this  thesis,  the  value 

of  v/c  is  expected  to  be  predictable  to  a  few  parts  in  10 

-4  -4 

Values  of  a    from  10    to  5  x  10    were  used  in  the  estima- 
tion presented  here.   It  is  expected  that  future  analysis 
of  Omega  data  for  these  frequencies  with  the  idea  of  con- 
sidering predicted  phase  velocities  will  yield  better  values 
for  p  and  a.   It  is  not  unlikely  that  different  standard 
deviations  of  predicted  phase  velocities  will  result  for 
different  frequencies.   If  this  does  occur,  K   in  (50)  can 
easily  be  modified. 

With  the  assumptions  delineated  above,  the  estimation 
can  be  performed.   As  discussed  in  Chapter  IV,  the  estima- 
tion procedure  consists  of  assuming  a  priori  knowledge  of 
an  unambiguous,  nominally  360  nmi,  lane  which  contains  the 
correct  receiver  LOP.   This  information  must  be  provided  by 
an  ambiguity  resolution  technique  in  a  final  system  imple- 
mentation.  Then,  the  sum  of  the  great  circle  distances  from 
the  two  stations  to  some  point  in  the  region  of  the  receiver 
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is  fixed  equal  to  L.   Since  two  stations  can  yield  only  an 
LOP,  this  sum  of  distances  can  be  quite  approximate,  as  will 
be  the  case,  since  this  sum  must  actually  be  obtained  from  a 
coarse  ambiguity  resolution  technique.   Finally,  x  is  de- 
fined as  the  distance  from  a  reference  station  (Forestport 
was  chosen  for  this  processing)  to  a  potential  receiver  site. 
Evaluation  of  the  likelihood  function  as  x  is  varied  over 
the  known  unambiguous  lane  provides  a  function  which  has  its 
maximum  value  for  that  x  which  is  the  maximum  likelihood 
estimate.   This  estimate  of  x  (xml)  is  equivalent  to  estima- 
tion of  an  LOP  as  can  be  seen  from  Figure  6.   The  path  over 
which  the  evaluation  is  done  is  shown  by  the  dashed  line,  a 
convenient  path  for  the  following  reasons.   First,  only  one 
parameter,  x,  need  be  varied  in  the  estimation  procedure; 
and  second,  for  a  fixed  step  size  in  the  evaluation,  the 
total  range  over  which  x  is  varied  (one  unambiguous  baseline 
lane  width)  is  the  same  for  any  receiver  position.   This 
method  of  estimation  necessarily  yields  all  of  the  local 
maxima  of  the  likelihood  function.   Thus  a  selected  number 
of  estimates  of  x  which  are  less  likely  to  be  the  correct 

estimate  than  x  n  can  easily  be  retained.   As  discussed  in 

ml 

Chapter  IV,  the  ratio  of  the  likelihood  function  evaluated 
at  these  alternate  candidate  values  of  x  divided  by  the  like- 
lihood function  evaluated  at  x  .  is  a  measure  of  the  relative 
probability  of  the  respective  values  of  x  being  correct. 

Figure  13  shows  a  block  diagram  of  the  estimation  pro- 
cedure.  The  collected  Omega  data  is  provided  to  a  processor, 
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described  in  detail  in  Chapter  IV,  which  evaluates  a  likeli- 
hood function  over  a  prescribed  range  of  values  for  x.   Para- 
meters which  are  provided  to  the  processor  are  the  mean 
values  and  covariance  matrix  for  the  phase  velocities. 

The  IBM  360-67  at  the  Naval  Postgraduate  School  was 
utilized  for  the  processing  of  the  Omega  four-frequency  data. 
A  Fortran  program  (Appendix  E)  ,  requiring  34K  BYTES  of  core, 
was  used  to  iteratively  evaluate  the  likelihood  function  for 
x  at  0.5  nmi  increments  from  x  =  855  nmi  to  x  =  1540  nmi. 
The  one  nmi  regions  centered  at  the  values  of  x  which  yielded 
the  three   largest  values  for  the  likelihood  function  were 
then  iteratively  searched  with  0.05  nmi  increments  to  obtain 
a  fine  estimate.   The  values  of  the  likelihood  function  for 
the  second  and  third  most  likely  values  of  x  were  then 
divided  by  the  value  of  the  likelihood  function  for  x  ,, 
yielding  a  probability  ratio  discussed  earlier.   The  time 
required  to  process  one  run  (a  data  set  corresponding  to 
three  minutes  of  collected  data)  was  on  the  order  of  15 
seconds  of  computer  time.   The  range  of  x  for  which  the  pro- 
cessing was  done  was  nearly  twice  the  nominal  360  nmi  which 
would  have  been  required  had  an  ambiguity  resolution  tech- 
nique provided  the  unambiguous  lane  over  which  x  was  to  be 


Any  number  of  values  may  be  retained,  however,  three 
were  found  to  be  adequate  to  always  contain  the  correct 
distance . 
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varied.   This  was  considered  the  only  fair  way  to  process 
the  data  since  the  correct  value  of  x  (1197.53  nmi)  was 
known  a  priori  and  was  considered  as  the  center  of  the  un- 
ambiguous lane.   In  a  finally  implemented  system,  the  a 
priori  knowledge  of  x  will  consist  of  some  probability  den- 
sity of  x,  across  one  or  more  unambiguous  lane  widths. 

C.   ESTIMATION  USING  PREDICTED  FIRST  MODE  PHASE  VELOCITIES 

As  an  intermediate  step,  the  optimum  estimator  must  esti- 
mate the  amplitudes,  noise  density  heights  and  doppler  fre- 
quencies of  all  eight  signals.   Appendices  A  and  B  give  the 
derivation  of  these  estimators.   The  doppler  shift  of  the 
signals,  which  would  result  from  physical  motion  or  frequency 
drifts  in  the  receiver,  were  found  to  be  negligible  since 
there  was  no  SARCOM  and  thus  no  motion  and  were  assumed  to 
be  zero  for  this  experiment.   In  a  finally  implemented  system, 
doppler  effects  must  be  estimated  and  removed.   The  signal 
amplitudes  and  noise  density  heights  which  were  estimated  as 
an  integral  part  of  the  processing  were  used  to  form  the  sig- 
nal power-to-noise  density  ratios   in  dB-Hz.   The  results 
of  this  estimation  are  presented  in  Appendix  B. 

The  results  of  the  LOP  estimation  utilizing  the  optimum 
estimator  are  presented  in  Table  VI.   The  mean  values  of 
phase  velocity  shown  in  Table  II,  and  the  values  for  K   given 
by  (49)  and  (50)  were  used  in  the  processing.   Several  values 

of  a,  the  standard  deviation  of  the  phase  velocity  to  velocity 

-4 
of  light  (v/c)  ratios,  were  used,  but  only  a  =  2  x  10   ,  a  = 
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3  x  10   ,  and  a  =  4  x  10    are  shown.   The  results  of  using 

other  values  of  O    were  significantly  worse  than  for  those 
shown.   The  minor  lane  errors  which  are  referred  to  in  Table 
VI  are  defined  as  runs  for  which  one  of  the  peaks  of  the 
likelihood  function  adjacent  to  the  correct  peak  is  the  lar- 
gest value  of  the  function.   This  results  in  an  LOP  estima- 
tion error  which  is  nominally  6-7  nmi  from  the  correct  LOP 
on  the  baseline.   Major  lane  errors  would  be  defined  as  lane 
errors  greater  than  a  minor  lane  error.   However,  in  this 
experiment  using  the  values  of  a    stated  above,  major  lane 
errors  did  not  occur.   In  Table  VI,  the  column  labeled  "per- 
iod" contains  the  numbers  1  or  2  and  the  letters  D,  N,  or  T. 
These  refer  to  data  collection  periods  1  or  2  and  Day,  Night, 
or  Transition  condition  paths,  respectively. 

The  means  and  standard  deviations  of  the  LOP  estimation 
error  are  in  baseline  measure.   The  hyperbolas  which  are 
formed  by  the  isophase  lines  from  a  station  pair  in  a  phase 
comparison  system  diverge  off  the  baseline  (great  circle  pass- 
ing through  the  two  stations).   The  divergence  occurs  directly 
as  esc  (0/2),  the  "spreading  factor,"  where  9  is  the  angle 
formed  by  the  two  great  circles  passing  through  the  receiver 
location  and  the  stations  as  discussed  in  Chapter  I.   There- 
fore, phase  differences  at  points  off  the  baseline  correspond 
to  greater  distances  than  on  the  baseline.   The  spreading 
factor  for  the  Dallas  receiver  location  for  this  experiment 
was  1.88,  thus  the  distances  in  Table  VI  must  be  multiplied 
by  1.88  to  obtain  the  true  measure  for  Dallas.   The  practice 
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of  using  baseline  measure  was  selected  so  that  the  results 
can  be  more  readily  applied  to  the  general  case  where  a 
spreading  factor  of  from  1  to  perhaps  2  might  result. 

The  values  of  a  which  yielded  the  fewest  errors  were 

different  for  the  two  data  collection  periods.   The  value 

-4 
of  O    =  2  x  10    resulted  in  the  fewest  minor  lane  errors  for 

-4 
the  first  period  and  a  =  4  x  10    resulted  in  the  fewest  for 

the  second  period.   Increasing  the  value  of  a  is  equivalent 
to  desensitizing  the  estimator  to  the  mean  values  of  phase 
velocity  being  used.   It  would  thus  appear  that  for  the 
data  obtained  in  the  second  period,  the  mean  values  or  pre- 
dicted values  of  phase  velocity  which  were  used  were  not  as 
close  to  the  actual  phase  velocities  which  were  encountered. 
It  is  believed,  however,  that  the  data  was  significantly 
disturbed  by  higher  order  modes  of  propagation  for  the  night 
path  from  the  Forestport  station.   Therefore,  use  of  the 
first  mode  phase  velocity  values  given  in  Table  II  would  not 
be  expected  to  yield  error-free  estimation  results  for 
night  paths.   This  conclusion  was  reached  after  the  follow- 
ing considerations.   The  mean  values  of  the  phase  differences, 
A<{>,  for  all  four  frequencies,  f,  for  all  lighted  (day)  trans- 
mission paths  were  utilized  to  estimate  v/c  from  the  data 
by  using 


(v/c) 


f.     (L-2x) 
c(n.    +   AcJ>  .  /2tt) 


Ki<4         (51) 
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to  obtain  Figure  14.    In  (51)  L  is  the  (known)  sum  of  the 
distances  from  the  stations  to  the  receiver,  x  is  the  (known) 
distance  from  the  reference  station  to  the  receiver,  c  is 
the  velocity  of  light  in  vacuum,  and  n   is  an  integer  equal 
to  the  number  of  whole  cycles  of  phase  difference  at  the  re- 
ceiver.  Equation  (51)  assumes  that  the  phase  velocities  for 
the  same  frequency  signals  were  the  same  for  both  paths  under 
day  conditions.   This  is  considered  to  be  reasonable  since 
the  paths  were  both  predominantly  westerly  and  first  mode 
propagation  during  day  conditions  apparently  dominated.  This 
latter  assumption  seems  to  be  borne  out  by  the  fact  that  no 
day  condition  lane  errors  were  experienced.   The  apparent 
parallel  shift  of  the  curves  in  Figure  14  from  the  predicted 
first  mode  values  could  be  explained  by  timing  errors  be- 
tween the  two  stations.   Several  0.25  \is    and  0.5  ys  correc- 
tions in  timing  were  inserted  by  Forestport  during  the  data 
collection  periods.   This  parallel  shift  in  phase  velocities 
will  not  cause  lane  errors  to  result,  but  will  be  seen  as  an 
estimation  error  or  a  bias  in  the  estimation  error  mean  value 
for  this  experiment. 

The  results  of  Figure  14  can  now  be  used  to  obtain  phase 
velocity  values  for  dark  paths  (night)  for  all  four  frequen- 
cies from  each  station  separately  to  determine  if  the  phase 


The  predicted  phase  velocity  values  shown  in  Figures  14, 
15,  16,  and  17  are  taken  from  [11,  12,  and  16]  and  private 
communication  with  Professor  J.  A.  Pierce  of  Harvard  Univer- 
sity.  The  lines  drawn  between  the  plotted  points  serve  only 
as  aids  for  identification  of  points  and  should  not  be  con- 
strued as  portions  of  continuous  curves. 
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velocities  from  one  station  appear  to  be  causing  the  lane 
errors.   This  was  done  by  using  the  change  in  the  mean  values 
of  the  phase  differences  from  day  to  night  for  each  frequency 
from  each  station  separately  in  the  equation 


(v/c) 


2TTf.d   (v/c)±tD 


i,N 


C<*i,N  -  *i,D)(v/c)i,D  +  2*fid 


Ki<4  .   (52) 


In  (52)  d  is  the  distance  from  the  station  under  considera- 
tion to  the  receiver,  (v/c).    are  values  taken  from  Figure 

1 ,  u 

14,  and  (<b  .  „  -  <t>  .  ^)  is  the  difference  of  the  mean  values 
'        i, N     i,D 

of  the  night  phase  and  the  day  phase  for  the  ith  frequency, 
treating  the  stations  separately.  The  results  of  the  night 
phase  velocities  obtained  in  this  way  are  presented  in  Fig- 
ures 15  and  16  for  the  first  and  second  data  collection 

3 
periods  respectively.    The  phase  velocities  from  the  Forest- 
port  station  are,  by  this  analysis,  apparently  the  cause  of 
the  lane  errors  which  were  experienced  for  the  night  condi- 
tions.  The  relatively  parallel  shift  of  the  Trinidad  phase 
velocities  from  the  predicted  first  mode  values  indicates 
that  the  Trinidad  data  is  probably  dominated  by  first  mode 
propagation.   However,  the  erratic   shifting  of  the  Forest- 
port  phase  velocity  values  from  the  first  mode  predictions 
indicates  very  probably  that  higher  order  propagation  mode 
interference  was  experienced  during  the  night  conditions  for 
both  data  collection  periods. 


Ibid. 
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One  further  processing  of  the  data  was  performed  to 
demonstrate  the  consistency  of  the  estimator.   Phase  velocity 
values  for  each  frequency  for  lighted  transmission  paths  and 
for  each  frequency  for  dark  transmission  paths  were  deter- 
mined for  each  data  set  using  (51).   This  method  determines 
effective  phase  velocities,  assumed  the  same  from  each  sta- 
tion for  the  same  frequency  and  path  conditions.   The  phase 

velocities,  normalized  by  c,  thus  determined  are  shown  in 

3 
Figure  14  for  lighted  paths  and  Figure  17  for  dark  paths. 

The  median  values  of  phase  velocities  of  the  two  collection 
periods  were  then  taken  from  the  curves  and  the  entire  data 
set  was  processed  using  this  one  set  of  phase  velocity 
values . 

The  results  of  this  processing  were:  three  minor  lane 
errors  for  the  second  period  with  night  conditions  and  one 
minor  lane  error  for  the  second  period  with  transition  condi- 
tions.  All  other  data  were  lane  error-free.   This  demon- 
strates that  if  reasonably  accurate  effective  phase  velocity 
values  are  known,  extremely  good  estimation  results  can  be 
expected  from  the  maximum  likelihood  estimator.   It  further 
demonstrates  that  at  least  for  the  experiment  conducted,  the 
mean  values  of  the  effective  phase  velocities  did  not  change 
significantly  in  more  than  one  month.   A  single  phase  veloc- 
ity value  for  each  frequency  for  night  and  likewise  for  day 
conditions  yielded  excellent  results  for  data  collection 
periods  separated  by  a  month. 


Ibid. 
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As  stated  before,  higher  order  mode  interference  is  the 
suspected  cause  of  the  lane  errors  experienced  in  processing 
the  data  utilizing  the  phase  velocity  values  from  Table  II. 
This  and  the  demonstrated  relative  insensit ivi ty  of  the  esti- 
mator to  phase  velocity  values,  along  with  the  lane  error- 
free  day  path  condition  results,  indicate  that  very  good 
estimation  results  can  probably  be  expected  from  using  single 
phase  velocity  values  for  a  given  receiver  location.   The 
sensitivity  of  these  phase  velocity  values  to  receiver  loca- 
tion remains  to  be  determined. 

D.   ESTIMATION  USING  PROPAGATION  CORRECTION  TABLES 

The  Defense  Mapping  Agency  publishes  Omega  Propagation 
Correction  tables  regularly  for  the  active  Omega  stations  in 
pairs  for  the  two  frequencies  10.2  kHz  and  13.6  kHz.   These 
corrections,  when  added  to  the  observed  phase  difference  of 
the  appropriate  frequency  from  a  station  pair,  allow  a  map- 
ping to  be  performed  from  the  observed  phase  difference  to 
an  Omega  navigational  chart.   These  charts  are  plotted  for 
an  assumed  value  of  (v/c)  =  1.00261  for  all  frequencies. 
Many  factors  in  addition  to  diurnal  effects  which  affect 
phase,  or  equivalently  the  effective  phase  velocity,  of  the 
signals  is  accounted  for  by  the  computer  program  which  gen- 
erates these  tables.   In  order  to  determine  whether  the  pre- 
dicted corrections  would  improve  the  results  of  the  experiment, 
tables  for  the  appropriate  location,  time,  and  station  pair 
were  obtained  from  Defense  Mapping  Agency.   The  corrections 
used  for  the  10.88  kHz  and  11  1/3  kHz  signals  were  obtained 
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by  linear  interpolation  from  the  published  corrections. 
Curves  showing  24-hour  plots  of  these  correction  factors 
are  shown  in  Figures  18  and  19  for  the  second  half  of  Octo- 
ber 1973  and  in  Figures  20  and  21  for  the  second  half  of 
November  1973. 

The  results  of  the  estimation  when  the  data  were  cor- 
rected using  the  factors  found  as  described  above  and  using 
(v/c)  =  1.00261,  for  all  frequencies  is  presented  in  Table 
VII.   The  number  of  lane  errors  was  decreased  by  utilizing 
the  propagation  correction  factors;  however,  quite  a  large 
number  of  lane  errors  still  occurred  at  night.   It  is  be- 
lieved that  this  helps  to  substantiate  the  higher  order  mode 
interference  theory  of  the  last  section. 

E.   EXPERIMENT  ESTIMATION  SUMMARY   • 

Because  of  the  suspected  higher  order  mode  interference 
at  night  during  this  experiment,  it  cannot  be  considered 
conclusive  that  use  of  the  Table  II  phase  velocity  values 
are  inadequate  for  LOP  estimation  in  a  final  SAR  system. 
Further  experiments  at  points  an  adequate  distance  from  the 
stations  to  preclude  higher  order  mode  interference  may  show 
that  propagation  correction  factors  are  helpful  in  reducing 
lane  errors.   If  this  results,  the  corrections  can  quite 
easily  be  implemented  at  the  SARCEN  and  are  not  considered 
as  a  major  hindrance  in  implementing  a  successful  system. 

One  of  the  most  important  beneficial  factors  in  using 
the  maximum  likelihood  estimation  procedure  is  the  avail- 
ability of  multiple  estimates  along  with  their  associated 
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relative  probabilities  of  being  correct  which  are  obtained. 
When  lane  errors  were  made  using  the  phase  velocity  values 
from  Table  II  for  the  first  period,  the  correct  LOP  was 
always  the  second  most  likely.   For  the  second  period,  the 
correct  LOP  was  second  most  likely  for  all  but  eight  runs 
where  lane  errors  occurred.   For  seven  of  the  eight  runs  the 
correct  LOP  was  third  most  likely  and  for  one  run  it  was 
fourth  most  likely.   When  propagation  correction  tables  were 
used,  the  correct  LOP  was  second  most  likely  for  all  but  one 
fun  where  lane  errors  resulted.   For  this  run  the  correct 
LOP  was  third  most  likely. 

For  the  estimation  performed  with  the  Table  II  phase 

-4 
velocity  values  and  a  =  3  x  10   ,  the  probability  ratios  dis- 
cussed in  Chapter  IV  were  examined.   As  explained  in  Chapter 
IV,  the  ratios  of  the  likelihood  function  evaluated  for  two 
values  of  x  is  equivalent  to  a  ratio  of  the  probabilities  of 
these  x  being  the  correct  x  for  the  receiver  location.   Fig- 
ures 22-25  give  selected  probability  ratios  for  the  first 
data  collection  period  and  Figures  26-31  give  selected  pro- 
bability ratios  for  the  second  period.   The  symbol  ®  denotes 
the  probability  ratio  of  the  second  most  likely  x  when  the 
correct  x  was  the  most  likely.   The  symbol  A  denotes  the 
probability  ratio  of  the  second  most  likely  x  when  the  cor- 
rect x  was  second  most  likely.   The  symbol  | ♦ |  denotes  the 
probability  ratio  of  the  correct  x  when  it  was  third  or 
fourth  most  likely,  the  latter  case  occurring  for  only  one 
run  in  the  entire  data  set. 
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From  these  results,  it  appears  that  a  threshold  could  be 
established  such  that  when  a  probability  ratio  of  an  x  was 
above  a  threshold,  it  would  be  considered  as  likely  and  SAR 
forces  could  be  sent  to  the  site.   It  is  not  suggested  that 
the  probability  ratio  values  be  treated  as  actually  having 
meaningful  numerical  value.   It  appears  that  the  maximum 
likelihood  estimator  places  quite  a  lot  more  confidence  in 
the  data  for  this  experiment  than  the  results  of  the  experi- 
ment would  suggest  is  reasonable.   Whether  these  values  will 
be  more  meaningful  as  numerical  values  when  receiver  sites 
where  higher  order  mode  propagation  is  not  a  factor  must  be 
determined  by  future  experiment. 

The  information  of  "ranked"  position  estimates  which  can 
be  provided  is  considered  to  be  of  great  importance  for  a 
SAR  system.   If  the  most  likely  position  estimate  does  not 
yield  a  recovery,  rather  than  spending  hours  or  even  days  in 
an  ever-widening  search,  the  SAR  forces  can  be  vectored  to 
the  next  most  likely  position. 


118 


t   \ 


VII.   MULTIPLE  STATION  POSITION  ESTIMATION 

A.   AN  OPTIMUM  POSITION  ESTIMATOR 

The  mathematical  model  to  be  used  for  the  multiple  sta- 
tion problem  will  be  the  same  as  that  for  the  two  station 
problem  presented  in  Chapter  II  except  that  the  vectors  will 
be  of  length  PK  instead  of  2K  where  P  is  the  number  of  Omega 
stations  from  which  data  is  received.   Arbitrarily  ordering 
the  stations  as  one  through  P,  the  first  K  components  of  the 
received  column  vector 

Y(t)  =  S(t)  +  N(t)  (53) 

will  be  from  station  one,  the  second  K  from  station  two,  etc. 
Here  again,  K  is  the  number  of  different  frequencies  common 
to  all  Omega  stations  for  use  in  hyperbolic  navigation.   In 
(53),  S(t)  and  N(t)  will  be  PK  dimensional  column  vectors  as 
well.   The  components  of  S(t)  will  represent  the  signal  por- 
tion of  the  received  vector  Y(t)  and  the  components  of  N(t) 
will  represent  sample  functions  from  white,  Gaussian  noise 
processes,  the  ith  component  n.(t)  having  two-sided  noise 
density  height  N./2.   The  ith  signal  component  can  be  written 


as 


s.(t)  =  A.  cos  [(to.  +  A.)t  -  d>.  +  6.]  (54) 

X  1  1       1        rl       1 
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as  in  Chapter  II.   The  same  arguments  as  followed  In  Chapter 
II  allow  the  following  observations 


to  . 


0) 


K+i 


U 


2K+i 


0) 


(P-l)K+i 


(55) 


and 


9K+i  =  92K+i 


e 


(P-l)K+i 


(56) 


where  0.  (1  <^  i  <^  K)  are  independent  random  variables,  uni- 
formly distributed  on  (-IT,  TT  ]  . 

It  will  not  be  possible  to  follow  the  convention  in 
Chapter  II  for  defining  the  propagation  phase  delays.   There- 
fore, the  following  definition  will  be  made 


*pK+l  =  W 


•/ 


P+l 


dy 


pK+1 


7— r-  -  2fTm  _,,.  (Ki<K) 
(y)        pK+i    


(57) 


[o<p<(P-l)] 


where  v.  is  the  kth  phase  velocity  and  x   is  the  great  circle 

k  c  J  n  b 

distance  from  the  nth  Omega  station  to  the  receiver  location 
being  estimated. 

Similar  to  the  presupposed  unambiguous  lane  in  Chapter 
IV,  a  presupposed  area  on  the  globe  is  required  to  begin  the 
multiple  station  estimation  process.   This  must  be  obtained 
by  some  ambiguity  resolution  technique  using  the  data  re- 
ceived from  P  stations.   The  estimation  problem  now  consists 
of  evaluating  a  likelihood  function  for  all  points  on  the 
presupposed  area  and  determining  the  column  vector  x  which 
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maximizes  the  function  on  this  area.   The  likelihood  func- 
tion will  appear  as  a  rubber  sheet  spread  over  the  pre- 
supposed area  with  many  spikes  under  it  causing  many  "hills" 
in  the  sheet.   The  value  of  the  column  vector  x  which  cor- 
responds to  the  "highest  hill"  will  be  the  maximum  likeli- 
hood estimate  of  x  and  the  components  of  this  x  will 
uniquely  determine  the  great  circle  distances  from  the  Omega 
stations  to  the  maximum  likelihood  estimate  of  the  Omega 
receiver  position. 

Following  the  notation  and  arguments  of  Chapter  IV,  the 
likelihood  function  for  the  multiple  station  problem  can  be 
written  as 


/ 


K   P-l 


A(X)    Ev    f(6)  exP  {  I    Z       [R      cos(*      -  y     ) ] 

1=1  p=0 


cos  8   +   Z    E   [R      sin  (<fr      -  y     )]sin  9.}d9.  (58) 

1=1  p=0  r  r 


In  (58), 


R.2  =  (C.2  +  S.2)  (2A./N.)2   (1  <  i  <  PK) 
l       i      ill        —    — 


(59) 


and 


y.  =  tan"   (S./C.) 
*i  i   i 


(1  <  i  <  PK) 


(60) 


where  S.  and  C.  are  the  data  from  the  ith  signal  as  in  Chaptei 
IV.   Ev(*)  denotes  the  expected  value  with  respect  to  the  PK 
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dimensional  phase  velocity  vector  with  jointly  Gaussian  com- 
ponents as  described  in  Chapter  II.  Utilizing  Equation  (15) 
from  Chapter  IV,  some  algebra  and  some  trigonometric  identi- 
ties, (58)  can  be  written  as 


P-l 


A(X) 


11   I„  {[  Z   R^ 
P  = 


P-l  P-l 

+  2   Z    Z 


V     ,    •    *  -n   |.l-  !  ■       .,  RpK+i  RqK  +  i  COS 

x=l       p=0  p=0  q=p+l   r 


^i,P,q]1/2}    (q  -^^  ' 


(61) 


In  (61), 


til .      =  principal  value  of  (u  ,,  ,  .  -  u  „,  . 
ri,p,q    r  r  qK+i     pK+i 


vqK+i    TpK+iy 


(62) 


Arguments  similar  to  those  made  in  Chapter  IV  allow  the 
high  SNR   maximum  likelihood  estimation  equation  for  multiple 
stations  to  be  written  as 


K      P-l    P-l 
A(X)    =    E       exp[    ZEE  Q.  i|> .  ]     . 

£.  i  =  l    p  =  0    q  =  P+l       X'P»1       X'P'« 


(63) 


In  (63)  e  is  the  PK  dimensional  column  vector  with  jointly 
Gaussian  components  just  as  defined  in  Chapter  IV  and 


Q. 
i,p,q 


R    .  R  „  . 
pK-fi   qK+i 

2(RpK+i  +  RqK+l> 


(64) 
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The  matrix  form  of  Equation  (63)  could  be  formed  and  the 
required  expected  value  taken  as  in  Chapter  IV.   The  matrices 
involved  will  be  quite  large,  however,  since  in  general,  four, 
five,  or  even  more  Omega  station's  signals  are  expected  to  be 
received.   Since  the  solution  to  the  estimation  problem  in- 
volves matrix  inversions  for  each  set  of  x   (1  <  p  <  P)  ex- 

P     -    ~ 

amined   in  an  iterative  process,  the  processing  time  required 
for  the  implementation  of  the  optimum  position  estimator  will 
be  prohibitive.   Although  non-optimum,  the  most  reasonable 
multiple  station  position  estimation  processor  to  implement 
is  considered  to  be  a  pairwise  station  processor. 

B.   MULTIPLE  LOP  POSITION  ESTIMATION 

Multiple  LOP  estimation  will  require  the  estimation  of 
the  set  of  most  likely  LOP  for  the  (P-l)  different  pairs  of 
the  P  received  stations.   Consideration  must  be  given  to  the 
spreading  factor  (defined  in  Chapter  I)  when  choosing  the 
stations  as  pairs.   The  smallest  spreading  factors  which  can 
be  achieved  are  desirable. 

Frenkel  in  [23]  presents  a  maximum  likelihood  position 
estimation  procedure  which  finds  the  optimum  position  esti- 
mate when  given  multiple  LOP  estimates.   Although  presented 


4 
Although  some  of  the  possible  candidate  components  for 

X  can  be  discarded  before  the  matrix  inversion  is  required 

by  utilizing  the  "coarse  check"  procedure  described  in 

Appendix  E,  the  processing  time  for  the  optimum  multiple 

station  estimation  will  still  be  prohibitively  large  for  a 

SAR  application  as  shown  in  Appendix  E. 
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as  an  ambiguity  resolution  technique,  the  procedure  applies 
directly  to  the  multiple  station  position  estimation  problem 
when  the  stations  are  treated  pairwise  as  proposed.   The  co- 
variance  matrix  required  for  the  position  estimation  proce- 
dure can  be  calculated  as  in  [23]  if  the  variance  of  the 
phase  estimates  corresponding  to  the  LOP  estimates  are  known, 
In  order  to  approximate  what  these  variances  will  be,  we 
recall  Equations  (25a,  25b) 


0)  ,d 

a  Z    2 e   _  2iTm.         (1  <  i  <  K)     (65) 

rivicx       l  —    — 


where  d  is  the  great  circle  distance  of  the  receiver  from 

the  Omega  transmitter  generating  the  K  signals.   Since  e. 

(1  _<  i  _<  2K)  were  assumed  to  be  jointly  Gaussian  random 

variables  with  known  covariance,  the  variance  of  each  <t>  .  due 

i 

to  phase  velocity  uncertainty  only  is 

fui.d\2 


a.         =1 I   a 

*i     V  c  /     ei 


2  .  (66) 


It  has  been  shown  as  a  result  of  [17,  21]  that  the  noise  dis- 
turbances to  observed  phase  will  in  general  be  dominated  by 
the  phase  disturbances  due  to  phase  velocity  uncertainty. 
It  is  therefore  a  reasonable  assumption  to  consider  only 
phase  velocity  contributions  in  determining  an  idea  of  the 
variance  of  the  2K  phases  in  the  two  station  LOP  estimation 
problem. 
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As  a  result  of  [17]  and  the  experimental  results  re- 
ported on  in  Chapter  VI,  standard  deviations  of  the  phase 

-4 
velocity  ratios  (v/c)  on  the  order  of  3  x  10    are  expected 

for  the  four  frequencies  10.2  kHz,  11  1/3  kHz,  13.6  kHz, 
and  10.88  kHz.   Thus,  Table  VIII  can  be  created,  where  d  is 
the  great  circle  distance  from  the  station  generating  the 
signals  to  the  receiver  in  nmi.   If  the  standard  deviations 
of  the  four  frequencies  given  by  Table  VIII  and  the  corre- 
lation coefficients  given  by  Equation  (49)  in  Chapter  VI 
are  used  and  the  phases  are  considered  to  be  equally  weighted 
by  the  estimator,  it  can  be  shown  that 


CT./d  ~  .1253  x  10 


-3 


(67) 


The  value  of  a,  in  (67)  corresponds  to  the  standard  devia- 
tion of  the  phase  estimate  for  one  station.   The  standard 
deviation  corresponding  to  the  LOP  estimate  resulting  from 

a  station  pair  would  be  the  rms  of  the  a,  associated  with 

.  <t> 

the  two  stations.   The  value  of  d  required  in  (67)  can  be 
quite  approximate  and  can  be  the  great  circle  distance  in 
nmi  from  the  appropriate  Omega  station  to  the  center  of  the 
presupposed  area  on  the  globe  which  was  assumed  given  by  an 
ambiguity  resolution  technique.  The  value  of  a,     given  by 
(67)  is  approximate,  since  the  optimum  estimator  actually 
accounts  for  the  estimated  SNR  of  the  signals  instead  of 
weighting  them  equally.   However,  (67)  yields  a  result  which 
is  most  probably  adequate  to  enter  the  multiple  LOP  position 
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TABLE  VIII 


STANDARD  DEVIATION  OF  PHASE  DIFFERENCES 


i 

f . (kHz) 

d>i 

(d  in  nmi) 

1,5 

10.2 

.1188 

x  10 

2,6 

11  1/3 

.1320 

x  10~3 

3,7 

13.6 

.1584 

x  10"3 

4,8 

10.88 

.1267 

x  10~3 
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estimation  procedure  of  [23].   Note  that  this  o,  only  applies 
for  four-frequency  Omega  using  the  frequencies  specified  in 
Table  VIII. 

A  final  piece  of  information  which  must  be  provided  to 
the  multiple  LOP  position  estimation  procedure  in  [23]  in 
order  to  obtain  valid  results  is  an  effective  phase  velocity 
associated  with  the  phase  <j) .   This  is  required  in  order  to 
do  the  mapping  from  phase  to  lat itude/ longitude  on  the  earth. 
Although  the  phase  velocities  of  the  different  signals  have 
be'en  modeled  as  having  known  mean  values  with  additive  dis- 
turbances, it  is  expected  that  an  average  of  the  mean  values 
of  phase  velocity  will  be  sufficient  to  provide  the  necessary 
phase  velocity  for  the  mapping.   These  mean  values  will 
already  account  for  diurnal  effects  on  the  path  being  con- 
sidered . 
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VIII.   CONCLUSIONS 

This  dissertation  has  reported  the  results  of  research 
efforts  applied  to  the  problem  of  achieving  an  optimum  posi- 
tion estimator  to  be  utilized  for  the  GRAN  worldwide  SAR 
system.   The  principal  new  results  which  it  contains  are 
discussed  below.   Areas  into  which  future  research  is  needed 
are  discussed  as  well. 

A.   PRINCIPAL  NEW  RESULTS 

The  dissertation  has  presented  a  new  way  to  model  VLF 
phase  uncertainties.   As  presented  in  Chapter  II,  these  uncer- 
tainties have  been  modeled  as  effective  phase  velocities 
which  are  jointly  Gaussian  random  variables.   The  propagation 
effects  on  phase  which  can  be  reasonably  well  predicted  such 
as  diurnal  effects  are  accounted  for  by  specifying  the  mean 
values  of  the  phase  velocities.   All  the  remaining  effects 
which  cannot  be  so  readily  predicted  are  modeled  as  contribut- 
ing to  the  uncertainty  in  phase  through  a  randomness  of  effec- 
tive phase  velocity.   This  model  allows  for  very  conveniently 
incorporating  any  correlations  in  the  phase  uncertainties, 
either  in ter frequency  or  inter station .   This  is  accomplished 
by  means  of  treating  the  effective  phase  velocities  as  being 
jointly  Gaussian  with  known  covariance. 

The  maximum  likelihood  LOP  estimator  for  an  arbitrary 
number  of  Omega  frequencies  was  derived  for  two  stations  in 
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Chapter  IV.   This  estimator  takes  advantage  of  the  model 
previously  discussed.   In  particular,  the  correlation  between 
frequencies  which  is  known  to  be  significant  and  is  presently 
actively  being  researched  is  simply  incorporated  into  the 
estimator.   If  the  inter frequency  correlation  which  is  ex- 
pected can  be  determined  quite  accurately,  lane  error  occur- 
rence can  be  very  significantly  reduced.   Uncertainties  in 
phase  or  equivalently  in  phase  velocity  would  still  cause 
errors  in  LOP  estimation  but  much  fewer  lane  errors  would 
occur  if  correlation  between  frequencies  from  an  Omega  sta- 
tion were  high  and  predictable. 

As  an  integral  part  of  the  maximum  likelihood  estimator, 
doppler  shift,  signal  amplitude,  and  noise  density  height 
estimation  equations  are  derived  in  the  appendices.   The  re- 
sults of  four-frequency  experimental  data  processed  by  the 
signal  amplitude  and  noise  density  estimators  are  presented 
in  terms  of  SNR  in  one  Hz  in  Appendix  B. 

Monte  Carlo  simulations  were  carried  out  utilizing  the 
maximum  likelihood  estimator  derived  in  Chapter  IV.   These 
simulations  were  carried  out  for  four-frequency  Omega  with 
fourth  frequencies  chosen  primarily  on  the  basis  of  ease  of 
implementation  at  the  Omega  stations,  their  separation  from 
other  known  or  planned  frequencies  and  the  resulting  unambi- 
guous lane  size.   As  discussed  in  the  previous  chapters,  it 
is  nearly  certain  that  for  any  of  the  known  methods  of 
ambiguity  resolution  which  appear  to  have  any  hope  of 
materializing  for  SAR  systems,  unambiguous  lane  size  in 
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excess  of  the  three-frequency,  nominally  72  nmi  lanes  must 
be  achieved.   This  fact  prompted  the  simulations  and  it  is 
now  possible  for  the  decision  concerning  the  choice  of  a 
fourth  frequency  to  be  made  in  terms  of  lane  error  tradeoffs 
to  achieve  acceptable  lane  size.   The  fact  that  the  simula- 
tions present  results  illustrating  the  extent  to  which  lane 
errors  can  be  expected  as  a  function  of  SNR  as  well  as  phase 
velocity  uncertainty  and  fourth  frequency  choice  allow  for 
observing  the  effect  of  SNR  on  lane  errors. 

The  simulation  results  prompt  the  conclusion  that  with 
adequate  knowledge  of  the  effective  phase  velocity  variance 
and  correlation  and  with  reasonable  signal-to-noise  ratios, 
the  unambiguous  lane  size  can  be  increased  beyond  the  present 
baseline  nominal  72  nmi  obtained  with  three  frequency  Omega. 
The  controlling  factors  in  the  selection  of  a  fourth  frequency 
for  Omega  to  be  used  for  SAR  purposes  are  (1)  the  capability 
of  ambiguity  resolution  techniques  to  resolve  the  lane  ambi- 
guity with  restricted  amounts  of  Omega  data  and  (2)  the 
accuracy  with  which  the  phase  velocities  of  the  Omege  fre- 
quencies can  be  predicted. 

The  maximum  likelihood  estimator  of  a  line-of-position 
for  the  two-station  Omega  problem  as  derived  in  Chapter  IV 
was  applied  to  experimental  four-frequency  Omega  data.   The 
results  of  this  estimation  (presented  in  Chapter  VI)  demon- 
strated that  the  maximum  likelihood  estimator  was  superior 
in  terms  of  lane  errors  to  the  walk-up  technique,  [29], 
which  is  presently  the  only  other  method  of  processing  the 
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experimental  data.   Of  particular  significance  was  the  20-25 
percent  fewer  lane  errors  experienced  using  the  maximum  like- 
lihood estimator  for  the  second  data  period,  night  path  con- 
ditions, over  the  walk-up  technique  employed  in  [29], 

The  experimental  data  processing  using  the  maximum  like- 
lihood estimator  was  done  utilizing  the  best  available  ex- 
pected inter frequency  effective  phase  velocity  correlation 
coefficients  (see  Chapter  VI)  and  various  values  of  standard 
deviation  for  these  phase  velocities.   The  standard  deviation 

values  which  yielded  the  best  results  in  terms  of  minimizing 

-4  -4 

lane  errors  were  on  the  order  of  2  x  10    to  4  x  10    as  ex- 
pected from  previous  study  by  other  authors  [17].   The  rela- 
tively large  number  of  lane  errors  which  were  experienced 
for  night  and  transition  periods  are  suspected  to  be  the 
result  of  higher-order-mode  propagation  interference  on  the 
Forestport  Omega  station  path.   The  basis  for  this  argument 
is  presented  in  Chapter  VI.   A  four-frequency  experiment 
planned  by  the  Naval  Air  Test  Center,  Patuxent  River,  Mary- 
land for  late  1974  which  is  to  relay  the  Omega  data  via 
satellite  from  various  globe  locations  is  expected  to  verify 
this  argument.   Retransmitter  (SARCOM)  positions  far  enough 
from  the  Omega  stations  to  insure  small  contribution  from 
higher-order-modes  of  propagation  are  planned. 

B.   FUTURE  AREAS  FOR  INVESTIGATION 

Several  areas  into  which  future  research  efforts  could 
be  made  have  been  observed  as  a  result  of  this  work.   Since 
this  is  the  first  known  attempt  at  modeling  the  VLF  phase 

131  i   \ 


velocities  as  random  variables,  prior  efforts  by  other  re- 
searchers have  not  analyzed  the  VLF  data  they  have  collected, 
in  a  fashion  which  readily  yields  the  necessary  parameters. 
The  variances  which  can  be  expected  as  a  function  of  dis- 
tance for  all  the  Omega  frequencies  need  to  be  better  deter- 
mined . 

There  is  presently  some  interest  in  investigating  the 
correlation  between  phase  uncertainties.   Professor  Pierce 
of  Harvard  has  been  examining  the  possibility  of  using  the 
inter  frequency  correlation  of  Omega  to  advantage  for  some 
time.   It  is  felt  by  this  author  that  the  determination  of 
correlation,  both  int er frequency  and  inters  tat  ion ,  if  it  can 
be  determined  to  be  significant  even  for  some  paths,  is  a 
most  fruitful  area  for  study.   Int erf r equency  correlation 
could  reduce  lane  error  occurrence  in  Omega  dramatically. 

In  terms  of  the  GRAN  application,  one  of  the  areas  which 
appears  to  have  been  neglected  is  the  study  of  methods  to 
reduce  doppler  shift  of  the  signals  received  at  the  SARCEN. 
Proper  design  of  the  SARCOM  and  processing  at  the  SARCEN 
could  eliminate  some  of  the  sources  of  doppler  shift  and  thus 
allow  for  an  estimate  of  the  velocity  of  the  SARCOM  to  be 
made.   In  the  present  advanced  OPLE  concept,  this  was  not 
considered  possible  [28],   For  a  SAR  system,  particularly  for 
processing  a  distress  incident  at  sea,  SARCOM  motion  during 
and  after  transmission  is  a  very  real  possibility.   It  will 
be  much  easier  to  locate  a  drifting  life  raft  minutes  or 
hours  after  notification  if  some  velocity  estimate  is  avail- 
able. 
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APPENDIX  A 

DOPPLER  SHIFT  ESTIMATION  FOR  THE  TWO-STATION 

K-FREQUENCY  CASE 

Ideally,  if  we  knew  the  frequency  of  the  received  Omega 

signals  at  the  SARCEN  precisely  (both  oo .  and  A.  for  all  i) , 

we  could  perform  the  necessary  operations  to  obtain  Equations 
(12a,  12b)  in  Chapter  IV.   They  are 

10(j-l)  +  T, 
C.  =  I  f  y.(t)  cos  [(oo.    A.)t]dt  (A-la) 


\±    -   Z    I  y±(t)  cos  [(03±  +  A±; 

j  =  1  Ao(j-l) 


and 

10(j-l)  +  T, 


M     r    VJ   '     i 

E     /      y±(t)  sin  [(w±  +  A.; 

j  =  1   Ao(j-l) 


S±  =   Z     /      y_.  (t)  sin  [(UK  +  A..)t]dt  .         (A-lb) 


3 


The  frequencies  to.  (1  _<  i  _<  K)  are  known  quite  well,  however, 
because  of  the  A/R  tone  which  is  created  in  the  SARCOM  and 
sent  with  the  Omega  data  to  the  SARCEN  and  because  the  frequen- 
cies at  the  Omega  stations  are  well  controlled.   Therefore, 
we  can  form 

10(j-l)  +  T, 

(A-2a) 


M  VVJ   '     i 

: .  '=   Z  /      y  .  (t)  cos  to.  tdt 

i    •  i  /       x  1 

3  =  1  *10(j-l) 


and 


M     """^  -+  Ti 
S±'=   Z     /      y,  (t)  sin  OKtdt  (A-2b) 


'=   Z     /      y  .  (t)  sin  to .  tdt 
j  =  1   •'lO(i-l) 


(j-D 
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at  the  SARCEN  from  the  received  data. 

In  order  to  form  the  likelihood  function  for  A 
(1  <_    i  <_  2K)  the  definitions 

10(j-l)  +  T. 


and 


10(j-l) 


y.(t)  cos  to  .  tdt 
'  i         1 


10(j-l)  +  T. 


10(j-l) 


y.(t)  sin  oi.tdt 


(A-3a) 


(A-3b) 


are  made.   Recalling  that 

y.(t)  =  A.  cos  [(to.  +  A.)t  -  <b.  +  9.]  +  n.(t),   (A-4) 
J  i        l         i     i      Ti     l      i    '   v 


the  probability  density  functions  of  C ! .  and  Si .  conditioned 

upon  knowing  9.  and  <J> .  can  be  found  as  follows, 
li 


10(j-l)  +  T± 

C!.     /      A.  cos  t(w.  +  A.  )  t-(j>.+0  .  ]  cos  co. 
ij    J  i       v  i     i    Yi   iJ      i 

10(j-l) 
10(j-l)  +  T, 


tdt 


10(j-l) 
can  be  written  as 


/n.(t)  cos  co. tdt 
1  1 


(A-5) 


A     10(j-l)  +  T± 


C.  =  -^    /      cos(A.t-cj).  +  6.)dt  +  N'  .  .     (A-6) 
ij     2    y  i   Ti     i7       cij 

10(j-l) 

where  the  second  order  frequency  term  is  disregarded  as  not 

being  in  the  passband  of  the  system  and 
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.10(j-l)  +  T± 


cij 


r  XJ  ± 

I  n  ( t)  cos  U) .  tdt 

J   10(j-l) 


(A-7) 


The  integration  in  Equation  (A-6)  can  be  carried  out  yield- 
ing 

C«    =  ^   (sin  [10  A.(j-1)  +  A.T.  -  *  +  6.] 
J       1 

-  sin  [10  A.(j-l)-  <J).  +  6.]}  +  N'  .  (A-8) 
i  J      ri     i       ci] 

which  can  also  be  written  as 


AT 

C!  .   =  — h. —  a.  cos  [v.  .  +  n.  ]  +  N'..  . 
xj       2     i        ij     'iJ     cij 


(A-9) 


In  (A-9) 


a . 

1 


sin  (A.T./2) 
11 

A.T./2 
l  1 


(A-10) 


ni 


-  d>.  +  0  . 
1     1 


(A-ll) 


and 


Y 


iJ 


10  A.(j-l)  + 


A.T. 
1  1 


(A-12) 


Recalling  from  Chapter  II  that  n.(t)  is  a  sample  function 

from  a  white  Gaussian  random  process  with  two-sided  spectral 

density  N./2,  and  since  N1 . .  is  formed  as  a  linear  function 
1                 cij 

of  n. (t) ,  N1 . .  will  be  a  Gaussian  random  variable.   N1 . .  will 
1       cij  cij 

have  mean  value  zero  and  variance 


10(j-l)  +  T. 


E[N  .  .  ] 
cij 


/ 


Ni      1 

■z —  6(t)-t  cos  w.idT 

2       2       i 


N.T. 
i  1 


10(j-l) 


(A-13) 
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where  the  second  order  frequency  term  is  again  disregarded. 
Therefore,  C!.,  conditioned  upon  knowing  <p       and  6  ,  is  a 
Gaussian  random  variable  with  mean  value 


AiTi 
E[C|j]  =  — —  a±'  cos  [y         +  n±] 


(A-1A) 


and  variance 


VAR  t  C !  .  ] 


i  1 


(A-15) 


It  can  similarly  be  shown  that  S! .  is  a  Gaussian  random 
variable  with  mean  value 


A.T. 

e[s!.]  = V1  sin  [y--  +  n.l 


(A-16) 


and  variance 


VAR  [  S!  .  ]  = 


N.T. 


(A-17) 


Using  the  conditional  probability  density  functions 
found  above  and  the  fact  that  Cl .  and  S ! .  can  readily  be 
shown  to  be  uncorrelated ,  and  thus,  independent  since  they 
are  Gaussian  random  variables,  the  joint  conditional  density 
function  of  C. .  and  S. .  for  the  ith  received  signal  can  be 


written    as 


M 


f (c! . ,s! . I  a,   e,  x,  v)  =     n     (  „  _  )    . 
J  J=l         1  1 


A    T 

exp   {"  W7T  [cij   "  "V1  ai  cos   (Yij  +  ni)]2 


i  i 
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2  AiTi 

[S'  +  -V1 


H1T±    ij 


a±    sin  (y    +  r\±)  ]     }  . 


(A-18) 


This  can  be  written  as 

M 


2k. a. 

x    x 


f(cj  ,s;   U.,  e.)  =    t     (^-)  exP{^i 

J  J  j  =1  x    i  i 


[C!  .     cos     (y  .  .    +    n  .  ) 


13 


ij  i 


-S 


Ij     s±"(V±j    +    n±>- 


A.T .a . 
i    i    x 


]-  -V(c!42   +   s!    2)} 


N .T.  v    ij 
x    x 


ij 


2       2 
MT.A.    a.  2A.a.  M 

/         2.M  ,  x    x       i     ,  r       i    x ,  V//-.I 

=     (W?TT)     eXp[ 2lT ]  exp{^— [cos  n         2    (C'j     cos    Y 

xx  i  x  j=l 

M 

S!.     sin    Y..)-    sin    n .        E     (S|.     cos    Y--     +    C!.     sin    Y--)]     _ 
ij  ij  i    j=1       iJ  iJ  ij  ij 


xj 


J 


M 

Z   (C 


!.2  +  s!.2)}  . 


N.T.  .  n    xj      Xj 
x  x  j=l     J       J 


(A-19) 


Now  utilizing  the  definition  of  the  likelihood  function 
as  given  in  Equation  (9)  of  Chapter  IV  and  again  following 
the  practice  of  not  redefining  the  likelihood  function  when 
a  monotonic  function  of  it  is  taken  because  of  the  equiva- 
lence of  the  likelihood  functions, 


A(V-/ 


MT.A.2a.2       2A.a.  M 

exp[ 1~ — ]  exp{  X    1[cos  n.  *  £  (C!   cos  y± 

x  i  j=l    J        J 

M 

!  .  sin  Y..)~  sin  n .   Z  (S|.  cos  Y..+  Cl.  sin  v .  . ) 1 }f  (n  .  ) dn . 
xj        ±3  x     v  xj       'ij    xj       'xj/J     'i'   'i 


(A-20) 
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In  (A-20)  f(n  )  refers  to  the  probability  density  function 
of  ri  .   The  random  variables  n   (1  <  1  <  2K)  are  uniformly 
distributed  on  the  interval  (-it,  it].   Therefore,  Equation 
(15)  in  Chapter  IV  can  be  utilized  to  write  (A-20)  as 

MT.A.2a  2       2A.a.F. 
A(A.)  exp  [-    12^    ]  Io  [   *  X  1]  (A-21) 


where 


M  M 

F.  =  [(  E   C!.  cos  v..  -   E   S! .  sin  v..) 
J=l  J=l  J 


M  M  2  1/2 

+  (  Z       S!   cos  Y±1  +   2   C»   sin  y       ) Z] X/Z.(A-22) 

j=l  j=l  J 


The  present  OPLE  ground  station  and  the  planned  SARCEN 
call  for  processing  the  received  signals  at  frequencies  which 
result  when  the  A/R  tone  is  used  for  demodulating  the  folded 
or  compacted  spectrum.   This  will  result  in  processing  fre- 
quencies of  less  than  1  kHz.   This  method  of  processing  re- 
sults in  doppler  shifts  due  to  the  local  oscillator  in  the 
SARCOM  which  will  not  be  removed.   Thus,  although  because  of 
motion  of  the  SARCOM  relative  to  the  stations,  the  elements 
of  the  first  K  and  last  K  sets  of  signals  will  be  highly  cor- 
related; this  correlation  is  not  taken  advantage  of  and  the 
doppler  shifts  are  computed  separately.   Additionally,  this 
method  of  processing,  because  of  the  reasonably  unstable  oscil- 
lator planned  for  the  SARCOM  units,  will  not  allow  velocity 
of  the  SARCOM  to  be  predicted  with  any  confidence  [28] „ 
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Reference  [28]  also  contains  a  doppler  estimation  procedure 

which  is  similar  to  the  results  of  this  appendix  although 

nonop t  imum. 

Evaluation  of  (A-21)  over  an  expected  range  of  doppler 

shifts  in  increments  for"  each  signal  will  yield  estimates 

of  the  doppler  shift  of  each  of  the  2K  signals.   These  values 

of  Aj  (1  <  i  <2K)  can  be  used  to  find  the  values  of  a. 
i   •  —   —  i 

(1  <  i  <  2K)  and  y..     (1  *<  i  <  2K,  1  <  j  £  M)  .   Substitutions 
of  these  values  into  the  equations 


C . .  =  C! .  cos  Y. . 


S  !  .  sin  Y •  •  » 


(A-23a) 


S..    =   —   C'..    sin    y..    +  —   S  !  .    cosy..     , 
ij         cXjl      iJ  iJ         ot±      ij  ij 


(A-23b) 


M 
C.  =   Z    C.  .  , 
1    j=l    ^ 


(A-24a) 


and 


M 
j  =  l 


ij 


(A-24b) 


for  all  i  (1  _<  i  _<  2K)  yields  the  desired  doppler  corrected 
data.   The  data  elements  C.  and  S.  can  be  shown,  following 
the  procedure  in  the  first  part  of  this  appendix,  to  be 
Gaussian  random  variables  with  expected  values 


e  (c±   |   e±,  (f^)  = 


and 


A.MT  . 

x2   x  cos  (-  <b±  +   e.) 


A.MT. 
E  (S.  |  6.,  (j).)  =-  -i^  sin  (-  <j>±  +  9±)  . 


(A-25a) 


(A-25b) 
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Their  variance  will  be 


N  MT  . 
VAR  (C   |  9  ,  cf>.)  =  VAR  (S±  |  8±,  $  )  =      *      (A-26) 

4a 

As  the  magnitude  of  the  doppler  shift  of  the  ith  signal  in- 

2 
creases  causing  A.T.  ■*  ±  2it,  a.   will  approach  zero  as  seen 

by  (A-10) .   The  variance  of  the  data  will  get  very  large  and 

sin  d 


d 


the  data  will  have  an  extremely  small  SNR.   A  standard 

curve  will  show  that  for  A.T.  =  ±  tt  ,  a .   =  4/tt2  *  0.4.   In 

11       '   l 

other  words,  since  T.  is  nominally  one  sec.  for  all  i,  dop- 
pler shifts  of  magnitude  <^  0.5  Hz  can  be  corrected  with  a 
deterioration  in  SNR  o  f  <^  4  dB.   A  doppler  shift  of  0.5  Hz 
corresponds  to  a  relative  velocity  of  in  excess  of  10,000 
nmi/hr  at  the  VLF  frequencies,  so  with  proper  processing  at 
the  ground  station  utilizing  the  A/R  tone,  doppler  shifts 
which  are  experienced  in  GRAN  should  be  able  to  be  estimated 
with  little  deterioration  in  SNR. 
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APPENDIX  B 
SIGNAL  AMPLITUDE  AND  NOISE  DENSITY  ESTIMATION 

The  maximum  likelihood  estimator  derived  in  Chapter  IV 
requires  an  estimate  of  the  signal  amplitude  and  noise  spec- 
tral density  height.   It  will  be  assumed  that  the  data  avail- 
able for  this  estimation  is  the  received  data  corrected  for 
doppler  shift  utilizing  the  estimation  results  of  Appendix  A. 
Therefore,  the  values  of  the  Gaussian  random  variables  C.  . 
and  S..  (1  £  i  £  2K,  1  <_   j  <_   M)  given  by  Equation  (A-23)  will 
be  the  data.   The  value  of  K  and  M,  as  before,  denote  the 
number  of  frequencies  and  received  Omega  pulses  respectively. 
The  random  variables  C..  and  S..  after  correction  for  doppler 
shift  will  have  conditional  mean  values  and  variance  given  by 


A.T. 
E(C..  !  6.,  (j).)  =  -V1  cos  (-<$>.    +    0.)  , 


(B-la) 


A.T. 
E(S±j  |  6i,  4>±)=-^j±    sin  (-cj)i  +  0±) 


(B-lb) 


and 


N.T. 

var(c.  .  |e .  ,<j>.)  =  var  (s. .  le .  ,$.)=  -i-7 

li '  i'Ti         N  ij '  i'Ti    .   2 

i 


(B-lc) 


In  (B-l),  A.  is  the  amplitude,  T.  the  pulse  duration,  cf>.  the 
propagation  phase  delay,  0.  the  random  phase  shift  and  N  ./2 
the  two-sided  noise  density  height  of  the  ith  signal 
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(1  <  i  <  2K).   Also  in  (B-l) 


a 


AT 
sin(-^-i) 

A..  T. 
i  1 


(B-2) 


where  A.  here  is  the  estimated  doppler  shift  of  the  ith  sig- 
nal.  The  estimate  of  the  value  of  A.  and  thus  a.  will  be 
considered  known  for  this  derivation  since  it  was  assumed 
as  the  first  processing  step. 

The  joint  conditional  probability  density  function  of  the 
data  from  the  ith  signal  is 


f(c     ,   s       |e.,   *.)  =   4^_-)M  . 


2a.    M 
exp  {-  — J~  [  E   (C 


A±T.  2 

,...  -  -j-   cos  n.)  + 
i  i  j=i    J 


M  A.T. 

E   (S.  .  + 


j  =  l 


ij 


-z —  sin  n.)  ]  i 


(B-3) 


where 


i     Ti     i 


(B-4) 


This  joint  conditional  density  function  can  be  written  as 


2a 


i   xM 


f(c..,  s..    |   e±,  *.)  =  C^)*  . 
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2         2    2 
2a  A   MT 

exp  {~  v~±  t?i  +  ~ A±  T±  c±  cos  n± 


+  A±  T   sin  n± ] } 


(B-5) 


where 


M       2       2 

K.   -  ^   (c.  .   +  s. .  )  , 
j=l     J 


M 


c4  =  z  c. . 

3  =  1 


(B-6a) 


(B-6b) 


and 


M 

£ 

3  =  1 


ij 


(B-6c) 


The  definition  of  the  likelihood  function  was  given  in 
Equation  (  9  )  of  Chapter  IV.   Again,  the  practice  of  not  re- 
defining the  likelihood  function  when  a  monotonic  function 
of  it  is  taken  will  be  followed  because  of  the  equivalence 
of  the  likelihood  functions.   The  likelihood  function  for 
amplitude  and  noise  density  estimation  can  thus  be  written  as 


A<W  -^b/ 


2  2    2 

2a.         A.  MT. 

i  i 


2a.    A. 
exp[ 1 (C±    cos    T]±    ~    S±    sin    n±) ]  • 


£(r\±)    dn±    . 


(B-7) 
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f  ( ri  )  in  (B-7)  is  the  probability  density  function  of  T).  . 
In  Chapter  II  0   (1  _<  i  _<  K)  were  shown  to  be  uniformly  dis- 
tributed on  (-tt,  tt]  ,  so  all  the  r| .  will  also  have  this  den- 
sity function.   This  and  Equation  (15)  in  Chapter  IV  allow 
the  likelihood  equation  to  be  written  as 


A(A±  N±) 


2  2    2 

2a.  A   MT 

— m  exp   ["  nTtT  (?i  + r~~)] 

N,  i  i 


2a.  A. B. 
z         (       x      x    x) 

o   ^    N,    ' 


(B-8) 


where 


2      2  1/2 
3.  =>  (C.   +  S  .  )  ' 

X  1        1 


(B-9) 


I  (•)  in  (B-8)  is  the  zero  order  modified  Bessel  function  of 
o 

the  first  kind. 

For  estimation  of  A.,  the  likelihood  function  is  written 

l 

as 


2      2 
A.  MT.a. 

A(A±)  =  exp  ( X2N  X  1  ) 


2a.  A. 3. 

I   (   X   X  X) 

o  v   N.     ' 


(B-10) 


Taking  the  partial  derivative  of  (B-10)  with  respect  to  A. 
and  setting  the  result  equal  to  zero  yields 


2~ 

2a.  A. 3. 

I   (— — ) 

2B  .     1  l    N.    ; 


i    MT 


2^ 
2a.  A. 3. 
I   ( 1 L_i) 

l 


(B-ll) 
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where  A.  is  the  estimate  of  A.  and  I, (•)  is  the  first  order 
x  i       1 

modified  Bessel  function  of  the  first  kind. 

For  estimation  of  N , ,  the  likelihood  function  is  written 
as 


A(N±) 


N 


2  2  2 

2a/  A.    MT. 

~m   exp[_  nTtt  a±  +  — — )] 


i  i 


2a.    A. 3. 

i    ( * tHi) 


(B-12) 


Taking  the  partial  derivative  of  (B-12)  with  respect  to  N 
and  setting  the  result  equal  to  zero  yields 


2a.  A. 3. 

2         2   2  2      1  ( — -) 

2a.  £.    a.  A.  T.    2a .  A .  3 .   1     N.    ; 

^  X    1   ,     1    1    X         X    XX  X 

1  "   MT.         2      "      M        „   2.  Q 

x  2a .  A . 3 . 

I  (   X   X  X) 

o^    N.    ' 


(B-13) 


where  N  is  the  estimate  of  N.. 

x 

Equations  (B-ll)  and  (B-13)  do  not  have  the  values  to  be 
estimated  isolated  and  further,  they  must  be  jointly  esti- 
mated since  each  estimate  is  dependent  upon  the  result  of 
the  other.   A  convenient  method  of  jointly  solving  the  equa- 
tions is  to  isolate  the  ratio  I..  (•)/!   (*)  in  each  equation 
and  equate  them.   This  yields 


N 


2      ~  2      2 
2a.  C-    A.   T.a. 

X   ^X       X      XX 

-      2 


(B-14) 


145 


(   \ 


Substitution  of  this  value  of  N.  into  (B-ll)  yields 


) 


4A.3.MT. 

I   (   1  1 i 

1  v      ~  2   2 

23,        45. -MA.  T. 

i  l    li 

MTi       4A.P.MT. 

!   (   1  1 i ) 

o         A  2   2 
°   4£.-MA.ZT.Z 
^i    li 


(B-15) 


Equation  (B-15)  can  be  iteratively  solved  for  A.  now  and  us- 

/\ 
ing  the  resultant  value  of  A.  in  (B-14)  yields  an  estimated 

<*. 

value  for  N  .  . 

l 

In  order  to  discover  a  method  convenient  to  solve  (B-15) 
iteratively,  the  following  was  done.   For  relatively  small 
N.,  in  other  words  high  SNR,  from  (B-14) 


2 
i   i 

MT. 


2m    2 
A.  T  .a. 
l   li 


(B-16) 


Or 


h      - 


2   2 
MA.  T. 
l   l 


(B-17) 


which  implies  the  argument  of  the  Bessel  functions  in  (B-15) 
are  large.   For  large  arguments,  I  (•)  and  I  (•)  are  approxi- 
mately equal  [27];  therefore,  for  high  SNR 


.   _  23. 

Ai   '  MT 


(B-18) 
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If  this  value  is  formed  from  the  data  for  each  i  (1  <_    i  <_    2K) 
and  defined  as  an  initial  condition  for  an  iterative  process 


i,o    MT.  * 


(B-19) 


the  (N+l)st  iteration  value  for  A.  can  be  shown  to  be 

x 


i,N+l 


=  (A,  J 


1,0 


h  < 


K    < 


">  2 
2A-  » 


2   2 
M  T. 

1 


4£.,    £.2„   MT.2 

1 ljN i_ 

"•2     2   2 
2A.  „   M  T. 

i,N 1 


) 


A  2 
4?,  "  A./ 

» J 


MT 


(B-20) 


A  computer  program  was  written  to  do  the  iterative  pro- 
cessing in  (B-20).   By  substituting  the  resulting  estimate 
of  A.  into  (B-14)  an  estimate  of  N.  was  obtained.   The  flow 
chart  of  the  computer  program  is  contained  in  Appendix  E. 
From  these  estimates,  the  signal-to-noise  ratios  in  one  Hz 
which  the  Omega  signals  have  at  the  front  end  of  the  SARCOM 
receiver  can  be  estimated  from 


SNR  = 


A2 


2N 


(B-21) 


The  results  of  evaluating  (B-21)  with  the  four-frequency 
data  obtained  from  the  experiment  discussed  in  Chapter  VI 
are  presented  in  Figures  Bl  -  B24.   Because  of  the  signifi- 
cant diurnal  effects  on  phase  velocity,  as  discussed  in 
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Chapter  II,  the  data  has  been  divided  into  three  major  cate- 
gories. Both  transmission  paths  in  daylight  is  referred  to 
as  "Day"  condition,  both  paths  dark  is  referred  to  as  "Night" 
condition  and  partially  lighted  paths  is  referred  to  as 
"Transition"  condition.   The  SNR  estimation  results  of  the 
first  five-day  data  period  were  so  little  different  from 
the  second  five-day  period  that  the  two  periods  have  been 
combined.   The  power  radiated  by  the  Forestport  Omega  station 
for  the  data  collection  periods  was  varied  over  quite  a  wide 
range  as  weather  conditions  varied.   Radiated  powers  on  the 
order  of  100-200  watts  for  the  lower  three  frequencies  and 
200-400  watts  for  13.6  kHz  were  typical.   The  radiated  power 
from  the  Trinidad  station  was  on  the  order  of  1  kW  for  all 
frequencies.   The  finally  implemented  Omega  stations  are 
expected  to  radiate  10  kW  at  all  frequencies  consistently. 
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APPENDIX  C 
LOP  ESTIMATION  FOR  KNOWN  PHASE  VELOCITIES 

In  developing  the  model  for  the  estimation  problem  to 
be  solved,  the  first  approach  which  was  used  was  to  assume  the 
phase  velocities  were  well  known.   The  values  of  phase  veloc- 
ities which  were  intended  for  use  were  those  given  in  Table 
II,  the  first  mode  propagation  values  obtained  from  Professor 
J.'  A.  Pierce  of  Harvard  University.   The  derivation  of  the 
likelihood  function  for  this  model  is  carried  out  in  this 
appendix.   Simulations  with  several  possible  fourth  frequen- 
cies were  carried  out  and  they  are  presented  here  to  serve  as 
an  upper  bound  for  the  performance  which  can  be  expected  from 
the  maximum  likelihood  estimator  when  the  uncertainties  in 
the  phase  velocities  are  accounted  for  in  the  derivation. 

The  difficulty  with  assuming  well  known  phase  velocities 
is  that  as  the  signal-to-noise  ratio  gets  reasonably  large, 
the  estimator  places  great  confidence  in  the  data  received. 
Since  the  additive  noise  is  the  only  factor  disturbing  the 
data  in  the  known  phase  velocity  model,  as  SNR  increases, 
very  good  data  would  be  expected.   Literature  search  and  some 
preliminary  experimental  data  processing  revealed  that  the 
assumption  of  known  phase  velocities  is  not  justifiable  in 
a  model  and  thus  the  model  described  in  Chapter  II  was  created, 
Because  of  its  value  in  bounding  the  performance  of  the  esti- 
mator and  in  obtaining  a  feel  for  the  problem,  the  known 
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phase  velocity  estimator  is  derived  herein. 

As  was  pointed  out  in  Chapter  IV,  Equation  (16)  can  be 
used  as  a  starting  point  for  the  known  phase  velocity  esti- 
mator derivation  since  no  approximations  or  operations  in- 
volving the  phase  velocities'  randomness  have  been  performed. 
Therefore  we  have  as  a  likelihood  function 


A(x)  -   T   Io  {[R±2  +  RK+.2  +  2R.  RK+.  cos  ^..]1/2}  (C-l) 
i  =  l 


where  the  variables  in  (C-l)  are  defined  as  in  Chapter  IV. 
For  high  SNR,  the  same  arguments  as  applied  in  Chapter  IV 
can  be  applied.   The  likelihood  equation  can  thus  be  written 
as  in  Chapter  IV,  Equation  (30)  as 


K 
A(x)  =  exp  [  Z      Q   iK  ] 
i=l 


(C-2) 


since  no  application  of  the  randomness  of  the  phase  velocities 
was  made  to  that  point.   In  (C-2)  Q.  and  if),  are  defined  as  in 
Chapter  IV,  Equations  (17)  and  (31).   The  log-likelihood  func- 
tion, an  equivalent  likelihood  function  to  maximize,  is 
convenient  to  form  here,  particularly  for  computer  evaluation 
to  avoid  overflow  problems.   Therefore, 


A(x) 


K 
i  =  l 


q.  i>. 

^1  r  1 


(C-3) 


Monte  Carlo  simulation  for  several  fourth  frequencies  which 
could  easily  be  added  to  the  Omega  format  was  carried  out 
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utilizing  the  estimation  Equation  (C-3).   The  results  of  the 
simulation  in  terms  of  the  probability  of  lane  error  occur- 
rence are  presented  in  Figure  CI.   The  definition  of  lane 
error  here  is  as  defined  in  Chapters  IV  &  V. 

Two  general  conclusions  can  be  drawn  from  the  simulation 
assuming  known  phase  velocities.   First,  since  there  is  a 
degradation  in  probability  of  lane  error  on  the  order  of  .1 
to  .2  in  going  from  unambiguous  lane  size  of  a  nominal  500 
nmi  to  approximately  1000  nmi,  the  latter  size  can  serve  as 
a  'soft  upper  bound  in  selecting  fourth  frequencies  to  do 
simulations  when  jointly  Gaussian  distributed  phase  veloci- 
ties are  assumed.   Second,  a  fourth  frequency  of  12.364  kHz 
which  results  in  a  nominal  785  nmi  unambiguous  lane  appears 
to  be  worse  in  terms  of  lane  error  probability  than  either 
10.462  kHz  or  10.737  kHz  which  result  in  nominal  927  nmi  and 
1356  nmi  unambiguous  lane  sizes  respectively.   This  implies 
that  frequencies  chosen  from  the  lower  part  of  the  10-14  kHz 
spectrum  can  be  expected  to  yield  better  estimation  results 
than  those  chosen  from  the  upper  part  of  the  spectrum  in 
terms  of  lane  error  probability.   This  can  be  explained  by 
the  fact  that  the  ith  frequency  being  higher  yields  adjacent 
maxima  of  the  ith  component  of  the  likelihood  function  which 
are  closer  together  and  thus  the  chances  of  small  perturba- 
tions in  the  phase  of  the  ith  signal  caused  by  noise  will 
more  readily  cause  the  peak  of  the  likelihood  function  asso- 
ciated with  the  incorrect  LOP  to  be  larger  than  the  correct 
peak . 
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APPENDIX  D 
EXPECTED  VALUE  OF  A  SPECIAL  FUNCTION 

Chapter  IV  derived  the  likelihood  equation  given  by  Equa- 
tion (32)  as 


T       T 
A(x)  =  E   exp  [e  Ac  +  B  e  +  D]  . 


(D-l) 


In  (D-l)  E  (•)  denotes  the  expected  value  with  respect  to  £ 
where  £  is  a  column  vector  with  2K  elements  jointly  Gaussian 
distributed  with  zero  mean  values  and  covariance  matrix  K  . 

£ 

A  is  a  2K  square  matrix,  B  is  a  column  vector  with  2K  elements, 
and  D  is  a  scalar.   The  expected  value  in  (D-l)  can  be  written 
as 


A(x)    = 


exp (D) 


(27T)K|K     |1/2 


/....  / 


I  1r      Ttr      "I 

exp{-   -^le    K£      £ 


2K 


-  2£TA£  -2BT£]}  d£  d£2K 


(D-2) 


where  |*|  denotes  the  determinant.   The  argument  of  the  expo- 
nential in  the  integral  defined  as  G  can  be  written  as 

T   —  1         T         T      T    — 1  TT 

G  =  E  K     £  -  2e  Ae  -  2B  £  =  £  (K     -  2A)e-B  e-e  B,  (D-3) 

which  can  then  be  written  as 
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G  =  eT(K  1    -2k) c    -  UT(K  -1  -2A)e-eT(K  -1  -2A)U       (D-4) 

t.  t-  C- 


where 


T     T     -1       -1 
U   -  B   (K     -  2A) 

e 


(D-5a) 


and  thus 


U  »  (K  ~1    -  2A)  X    B  . 


(D-5b) 


Equation  (D-4)  can  now  be  written  as 

G  .=  (e  -  U)T  (K  -1  -  2A)  (e  -  U)-  UT(K  -1  -  2A)  U       (D-6) 

by  completing  the  square. 

These  manipulations  allow  (D-2)  to  be  written  as 


A(x) 


exp  [D  +  -|uT(K  ~1    -    2A)  U] 
(27T)K  I  K£  |  1/2 


/    /    exP{-  |  [(e-U)T(K£  X    -2A)(e-V)]}d£1    ...de2K  . 


2K 


(D-7) 


The  integral  is  over  all  values  of  the  jointly  Gaussian  random 
vector  e    and  thus  equals 


(2tt)K   I  (K  "1  -  2A)_1  | 


-1  i 1/2 
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The  desired  expected  value  is  then 


|  (K^1  -  2A)"1  |1/2 
A(x)  =  s r-p. exp  [D+  ±IT(K   X  -2A)U] 

Ik  |1/2  2    e 


.(D-8) 


Matrix  manipulation  of  this  equation  allows  the  final  simpler 
result 


A(x) 


exp  [D  +  BT  K  (I  -  2AK  )  1  B] 


(I  "  2AK£) | 


1/2 


(D-9) 


to  be  written. 
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APPENDIX  E 
COMPUTER  IMPLEMENTATION  OF  TWO-STATION  LOP  ESTIMATION 

In  order  to  solve  the  estimation  equation  given  by  Equa- 
tion (38)  in  Chapter  IV  for  the  two-station,  four-frequency 
experimental  data  discussed  in  Chapter  VI,  a  computer  pro- 
gram was  written.   A  flow  chart  of  that  program  is  contained 
in  this  appendix.   On  the  IBM  360-67,  the  program  occupied 
3'4K  BYTES  of  core.   The  time  required  to  process  one  "run" 
of  data  was  on  the  order  of  15  seconds.   A  run,  as  defined 
previously,  consisted  of  nominally  three  minutes  of  Omega 
data.   The  range  of  values  over  which  x  was  allowed  to  vary, 
as  discussed  in  Chapter  VI  was  twice  that  which  would  be  re- 
quired in  an  operational  system. 

A  "coarse"  processing  which  was  performed  to  save  com- 
puter time  was  the  following.   The  values  of  iJj.(1  <  i  <  4) 
as  given  by  Equation  (17)  in  Chapter  IV  were  examined  with 
the  value  of  x  being  tried  in  the  iterative  procedure  before 
entering  the  actual  matrix  manipulation  portion  of  the  pro- 
gram.  The  magnitude  of  each  value  of  ty .    was  required  to  be 
less  than  1.5  rad.  before  the  entire  estimation  process  was 
carried  out.   This  is  equivalent  to  saying  that  noise  and 
phase  velocity  uncertainties  will  not  disturb  the  phases  by 
more  than  +  1.5  rad.  for  the  Dallas,  Texas  receiver  location, 

The  value  which  is  used  for  this  "coarse"  processing  can 
be  determined  for  the  ith  frequency  by  examining  the  amount 
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of  phase  shift  which  two  standard  deviations  of  phase  veloc- 
ity disturbance  from  each  station  in  opposite  senses  will 
cause.   To  determine  the  effect  of  a  phase  velocity  distur- 
bance on  phase  the  equation 


*K+i  -  *i 


,L-x 

id  .  ( 

1   v 


K+i 


— )  -  2TTm. 
v  .        i 

1 


(E-l) 


is  written.   The  values  of  <j>  are  the  propagation  phase  delays 

for  the  ith  frequency  signals,  oj  .  is  the  ith  angular  frequency, 

y.  are  the  phase  velocities,  m.  is  the  number  of  whole  cycles 

in  the  phase  difference  and  (L-x)  and  x  are  defined  as  in 

Chapters  IV  and  VI.   To  see  the  maximum  effect  of  perturbing 

phase  velocities,  let  Av  be  added  to  vTr  ,  .  and  subtracted  from 
r  '  K+x 

v.  and  for  simplicity,  assume  v.  =  vT, ,  .  =  v.   Thus, 
i  v  j  ■>  l     K+i  ' 


*K+i    ~    *i    ~    A*i    =    Wi     (vJi?Tv- 


^-t-)    -    2TTm. 

v    -    Av  i 


to  , 


±   (■ 


L-x 


x 


1    + 


Av 


1    - 


7—)    -    2Trm. 
Av'  i 


(E-2) 


where  At})  .  is  the  perturbation  in  the  phase  difference  caused 
by  Av.   This  can  be  written  as 


0) 


<(>„,.  -  <f>.  -  A<j> . 
TK+i    ri     ri 


[(L-x)(l  -  £2L)-  x(l  +  ^)]-  2^m. 


W .                   to  .L   . 
—  (L-2x)  -  2TTm. —  (— ) 

V  1       V      V 


(E-3) 


since  ( )  will  be  small.  Or, 

v  ' 
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A<f>. 


w  .L 

1 


Av 
v 


w  .L 

x 


Av 
c 


(E-4) 


For  L  =  3547.15,  which  is  the  For es tpor t-Dallas-Tr inidad 

i A v i  —  4 

distance,  co.  =  2tt  •  13.6  kHz  and  |  —  |  =  8  x  10   ;  the  magni- 
1  c 

tude  of  the  maximum  expected  phase  perturbation  is  found  to 
be 


.  ,,  ~  2tt  •  13.6  •  3547.15  ;  8  «  10 
A*'  =  161875.2 


-1 


1 . 5  rad  . 


(E-5) 


Since  the  |  A<}>  |  for  lower  frequencies  will  be  smaller,  this 
was  considered  a  reasonable  upper  bound  to  use  for  all  fre- 
quencies of  the  four- frequency  experimental  data. 

The  amount  of  savings  in  terms  of  the  number  of  values 
of  x  which  must  be  considered  when  a  coarse  search  such  as 
this  is  used  can  be  determined  as  follows.   Consider  the  ith 
frequency  phase  difference  as  x  is  varied  through  a  lane. 
Only  one  half  the  values  of  x  will  meet  the  coarse  criterion 
if  \ty  .   |  is  required  to  be  less  than  -r-  .   Since  there  are 
four  frequencies,  only  1/16  of  the  possible  x  will  meet  the 
coarse  criterion.   As  discussed  in  Chapter  VI,  +  350  nmi 
(nominally)  were  used  in  processing  the  four-frequency  experi- 
mental data  with  increments  of  x  =  0.5  nmi.   Thus,  only  about 
90  values  of  x  vice  about  1400  had  to  be  processed  through 
the  entire  estimation  procedure. 
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For  a  generally  implemented  system  such  as  GRAN,  it  would 
not  be  unusual  to  have  values  for  L  =  7000  nmi  or  greater,  so 
it  is  not  expected  that  this  processing  advantage  can  be 
generally  implemented.   Because  of  the  significant  time 
savings  it  affords,  it  may  well  be  worth  using  a  variable 
value  which  is  determined  after  a  rough  estimate  of  L  is 
known.   In  the  final  GRAN  system,  if  10.88  kHz  is  added  as 
a  fourth  frequency,  360  nmi  unambiguous  lanes  will  result. 
If  Ax  of  0.5  nmi  are  used,  720  values  of  x  must  be  processed. 
If'no  coarse  estimate  procedure  is  used,  using  the  ~  15- 
second  processing  time  for  twice  the  number  of  values  of  x 
as  a  guide,  about  two  minutes  will  be  required  to  estimate 
an  LOP  for  a  station  pair  using  a  program  similar  to  the 
one  in  this  appendix. 

Preferable  methods  of  search  for  the  values  of  x  which 
maximize  the  likelihood  function  may  exist,  such  as  [30]. 
Methods  other  than  the  exhaustive  search  method  discussed 
in  previous  chapters  of  this  dissertation  and  implemented 
in  the  following  computer  program  were  not  considered.   The 
exhaustive  search  technique  was  chosen  because  of  its  sim- 
plicity and  ease  of  implementation. 
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c(l)  =  0 


S(I)  =  0 


p(D  =  o 


j  =  1 


Read     3 
(l,J)th 
Data 
VElementy 

i 

c(D=c(i} 

+  ec(l,J) 

IT 


~*gQ~ 


<D 
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s(i)=s(i) 

+  es(l,J) 


£(D=/3lW2 
(i.jV+ea'tl.J 


J  =  J  +  1 


SQ(I)=(S2{I) 

+  c2d; 


:(D) 


X 


R(D  = 
SQ(I)/9T(I) 


i. 


Wl)=   tan" 
(S(I)/C(I)) 


K!> 


^L   y 

I 

.,    \  > 

_J 

J   . 

1 

l 

V(I)=FR7-CL(I 
+(1-FR1)CD(I) 

V(I)=FR2»CL(I. 
+(1-FR2)CD(I) 

A 

U 

V 
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e».    1  =  1+1 


Call 
NEST     4 


R(D  = 

2R(I)-SQ(I) 


I  =  1 


/<I)    =yu(l+4) 
-yu(l)+SWC(l) 


, i 


.<$. 


Q(D  = 
-R(l)-R(l+4) 


2R(l)+2R(l+4) 


<V 


E»    1 


<D 
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XF-XS 


IND=0 


s 


J  =  1 


X=XS+(J-1)AX   O 


I  =  1 


W(l)=/a(l)-4i: 


I 


'(H)  _     X 
Ml+4)    ^J). 


1=1+1         «a- 


< 


e,       1=1+1 


,-iz\ 


J=  J+1  ■  <3 


<3 


W(l)=W(l)+2/T 


t        \ 
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1=1+1         <J- 


I  =  1 


I  =  1 


D  =  0 


D=D+W2(I)Q(I)  <sj 


Z=2Qmw(I)' 

MD/c) 


B(I)=-X-Z 


J 


B(l+4)  = 
(L-X)-Z 


77-Q(D(^(l) 


X 


A(I,I)=X2ZZ 


4 


A(i+4,I+4) 

(L-X)^.ZZ 


© 
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*  The  operations  in  these  boxes 
are  matrix  operations.  The  sym- 
bol | •  |  denotes  the  determinant, 
I  denotes  the  identity  matrix 
and  B  denotes  the  transpose  of 
the  vector  B. 


U3=£B  Kf« 
(I-2AK€ )~1B 


/6  «a- 


U1=(Jl-2AKj)« 


< 


U=U3+D-Ln(U1 ) 


XHOLD(K)=X 


J. 


HOLD(K)=U 


0 


111 


IND=1 


A X= AX/10 


K  =  0 


K=K+1 


JL 


XS=XHt>LD(K)- 
10AX 


XF=XHOLD(K)+ 
10AX 


XHOLD(K)=0 


HOLD (K)=-1 000 


E> 


N=20 


1 


<P 
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M  =  1 


HOLD(M)  = 

HOLD(M) 
H0LD(1) 


O- 


M=M+1 
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The  parameters  to  be  read  in  are  the  scalars  XS,  XF,  AX, 
L,  FR1,  FR2,  and  IDIM  and  the  matrix  K  . 

The  invariants  to  be  set  up  are  the  scalars  it,  RNG  and  c 
(vel  of  light  in  vacuum)  and  the  vectors  T,  a),  CL,  CD 
and  RNG. 

The  (I,J)th  data  element  consists  of   ec(I,J)  and  es(I,J) 
where 

10(J-1)+T(I) 
ec(I,J)  =  /  y(I)  COS[o)(I)t]dt  and 

10(J-1) 


10(J-1)+T(I) 
ec(I,J)  =  /  y(I)  SIN[uo(I)  t]dt 
10(J-1) 

where  y(I)  is  the  Ith  component  of  the  received  signal 
vector . 

Subroutine  NEST  must  be  provided  the  vectors  R,  3  and  T 
and  returns  in  R  the  estimates  of  the  signal  amplitude 
to  noise  density  ratios. 

The  vector  SWC  must  be  provided  as  propagation  correction 
factors  for  each  frequency  in  the  form  of  station  two 
minus  station  one  corrections.   When  propagation  correc- 
tion factors  are  used,  all  entries  in  the  mean  value 
vectors  V  must  be  set  equal  to  the  charting  phase  velocity 
which  is  used  in  their  prediction  (1.00261*c).   When  no 
propagation  correction  factors  are  used,  the  vector  SWC 
must  have  all  zero  entries. 

Subroutine  ORDER  must  be  given  the  values  of  the  scalars 
U,  X  and  IDIM  and  it  returns  the  requested  number  (IDIM) 
of  most  likely  estimates  of  X  in  the  vector  XHOLD.   The 
corresponding  likelihood  function  values  in  the  vector 
HOLD  are  returned  also.   The  components  of  these  vectors 
are  in  descending  order  of  likelihood. 

The  vector  HOLD  at  this  point  contains  the  probability 
ratios . 
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DEFINITIONS 

XS-  starting  value  for  X  in  the  iteration 

XF-  final  value  for  x  in  the  iteration 

AX-  size  of  incremental  steps  in  X  to  be  used 

L-  fixed  value  for  L 

K  -  8-by-8  covariance  matrix  of  £ 

CL-     4  element  vector  for  mean  values  of  phase  velocity 
in  lighted  path 

CD-     4  element  vector  for  mean  values  of  phase  velocity 
in  dark  path 

T-  8  element  vector  for  Omega  pulse  durations 

FR1-  fraction  of  path  1  in  light 

FR2-  fraction  of  path  2  in  light 

to-  4  element  vector  of  radian  frequencies 

y-  8  element  vector  of  mean  values  of  phase  velocity 

to  be  computed  by  the  program 

SWC-    4  element  vector  of  propagation  correction  factors 

in  form  of  station  two  minus  station  one  corrections 
for  each  of  four  frequencies  (to  be  set  equal  to 
zero  if  no  corrections  are  to  be  used) 

RNG-    a  predetermined  constant  in  radians  equal  to  the 
maximum  deviation  of  the  phase  of  any  signal  from 
the  phase  resulting  from  mean  values  of  phase 
velocity.   These  deviations  are  due  to  phase  velocity 
uncertainties.   (1.5  rad  was  used  in  the  experiment 
data  processing.) 

IDIM-   The  number  of  most  likely  estimates  of  X  to  be 
retained 

HOLD-   vector  of  length  IDIM  of  values  of  the  likelihood 
function  to  be  retained 

XH0LD-  vector  of  length  IDIM  of  values  of  X  corresponding 
to  the  values  of  HOLD 

y-     8  element  vector  of  phase  angles  computed  by  the 
program 
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SUBROUTINE 
NEST  ' 
Inputs 


Set  Up 
Constants 


K  =  1 


I  =  1 


X 


ZZ  = 

2<?(K) 


9R2(K)-T2(K) 


Z=  36/(ZZ-l) 


E=Z/3.75 


<D 


R1=1/E 


4l> 
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*  Q=C1(C2+C3-R1+C4'R2+C5-R3+C6-R4.+ 
C7'R5+C8'R6+C9'R7+C10'R8) 


R>R2/E 


R4^R3/E 


R5=R4/E 


R6=R5/E 


R7=R6/E 


R8=R7/E 


Q=. 


V 


R2=E 


R>R2«E' 


0 
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«*  Q=D1(D2-R1+D3-R2+D4-R3+D5'RA+ 
D6-R5+D7-R6+D8-R7) 


V 

R4=R3*E2 

V 

R5=R4*E2 

W 

R6=R5  -E2 

1 

' 

R7=R6-E2 

1 

r 

Q=           ** 

R(K)=R(K)/Q 


A(K)=R(K) 


-+>■     r(k)=i 


1=1+1 


-E> 


0 
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ZN 


A2(K)'T(K) 


SNR(K)   = 

A2(K) 

2-ZN(K) 


*«-  Ms 
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Subroutine  NEST  computes  the  joint  maximum  likelihood 
estimates  of  the  signal  amplitudes  and  the  one-sided 
noise  spectral  density  heights  from  the  equations 

Mrn   .  P(I)     A2(D  T(I) 

N(I)  -  lf7j) y^—         and 


.    .          I  [2A(I)-SQ(I)/N(I) 
(j\    -       ai^i;  .  —± 

y  1UJ    I  [2A(I)-SQ(I)/N(I) 


In  these  equations  the  vectors  are  as  defined  in  the 
subroutine  NEST  or  the  MAIN  estimation  program.   I   and 
I   are  modified  Bessel  functions  of  the  first  kind  of 
order  one  and  zero  respectively.   The  ratio  of  these 
Bessel  functions  is  computed  by  using  their  polynomial 
approximations  found  in  Handbook  of  Mathematical  Functions 
by  M.  Abramowitz  and  I.  A.  Stegun,  Dover  Publications. 

The  constants  required  for  this  subroutine  are: 


CI   = 

2.50662828 

Dl 

= 

.53333333 

C2   - 

.39894228 

D2 

= 

1.00000000 

C3   = 

.05316616 

D3 

= 

1.75781160 

C4   = 

.01118812 

D4 

= 

-1.02993610 

C5   = 

-.00161279 

D5 

= 

.90497756 

C6   = 

.01920037 

D6 

= 

-.84749981 

C7   = 

-.04017321 

D7 

= 

.80494379 

C8   = 

.04837180 

D8 

= 

-.76696533 

C9   = 

-.02678420 

CIO  = 

.00532070 

It  has  been  found  experimentally  that  10  iterations  are 
adequate  for  reasonable  convergence  of  the  iterative 
process  for  data  from  the  OPLE   ground  station  in  SNR 
environments  which  can  be  expected  from  Omega. 

The  vectors  A,  ZN,  SNR  and  R  written  out  here  are  the 
estimates  of  signal  amplitude,  two-sided  noise  density 
height,  signal  to  noise  density  ratio  and  signal  ampli- 
tude to  noise  density  ratio  respectively. 
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VSUBROUTINE  / 
ORDER 1   / 
i Inputs*  / 


I  =   1 


*  HOLD,XHOLD,X,U,IDIM 


-6»- 


P=  X-XHOLD(I) 


1=1+1 


E> 


j  =  i 


rs£. 


<D 


j=i+i 


HOLD(l)=U 


XHOLD(l)=X 
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H0LD(I-1) 
=  HOLD(I) 


XHQLD(I-I) 

=  XHOLD(I) 


HOLD(l)=U 


XHOLD(l)=X 


<3= 


1=1+1 


O 


Subroutine  ORDER  receives  the  scalars  X,  U  (likelihood  function  value 
for  X)  and  IDIM.  ORDER  places  U  and  X  into  the  vectors  HOLD  and  XHOLD 
of  dimension  IDIM  in  their  appropriate  place.  That  is,  the  largest 
value  of  U  received  thus  far  in  the  "run"  occupies  the  first  component 
of  HOLD,  the  second  largest  U  occupies  the  second  component,  etc. 
The  components  of  XHOLD  are  the  values  of  X  which  correspond  to  these 
values  of  U. 

This  operation  insures  that  only  one  value  of  U  and  X  will  be  retained 
for  each  peak  of  the  likelihood  function. 
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